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Introduction 


The  realization  of  structures  that  do  not  scatter  electromagnetic  field,  i.e.  structures  that 
appear  invisible  to  EM  waves,  is  a  concept  that  has  been  investigated  theoretically  since  the 
1960s  when  the  possibility  of  a  plane  wave  passing  without  distortions  through  a  structure 
with  anisotropic  filling  was  investigated  [1].  Hard  surfaces  have  been  used  in  the  design  of 
supporting  struts  of  reflector  antenna  feeds,  i.e.  it  was  experimentally  demonstrated  that  it  is 
possible  to  considerably  reduce  the  scattered  field  for  one  angle  of  incidence  and  for  both 
polarizations  of  the  incident  wave  [2],  In  [3] -[6]  it  was  shown  that  the  scattered  field  of  a 
two-layer  dielectric  ellipsoid  is  zero  for  certain  combinations  of  permittivities.  To  obtain  that 
at  least  one  of  the  layers  should  have  relative  permittivity  smaller  than  one  (i.e.  local  negative 
polarizability,  which  is  inherent  to  plasmonic  materials  and  plasma-like  ENG  metamaterials). 
Several  other  concepts  for  obtaining  invisible  scatterers  were  proposed,  like  minimum 
scattering  antennas  and  active  scatterers  ([7],  [8]). 

The  focus  has  recently  shifted  to  the  possibility  of  cloaking  objects  using  a  metamaterial 
cover  ([8]-[15]).  The  main  principle  of  metamaterial  cloaking  approach  arises  from  coordinate 
transformations  which  transform  the  volume  (in  particular,  cylindrical  one)  into  a  shell 
(cloak),  making  the  space  inside  the  shell  concealed.  The  required  cloak  material  is  fully 
anisotropic  and,  in  addition,  the  constitutive  parameters  (tensors  s  and  p)  are  functions  of 
radial  coordinate.  Such  cloak  design  is  not  feasible  in  practice,  however  the  simplified  cloak 
designs  were  proposed  in  [16]  and  [17]  which  claim  to  work  for  one  specific  polarization  of 
the  incident  plane  wave  (TMZ  and  TEZ  polarizations,  respectively).  These  specific  cloaks 
benefit  from  simplifications  of  constitutive  parameters,  since  only  radial  component  of 
permeability  [16]  or  permittivity  [17]  needs  to  have  radial  variation,  which  is  obtainable  by 
using  split-ring  resonators  or  metal  wires  in  dielectric  host,  respectively. 

Experimental  results  ([16],  [19]  and  [20])  reveal  that  the  measured  scattered  field  was 
reduced  compared  to  the  object  without  a  cloak  (i.e.  the  cloak  working  principle  was  proved), 
however  the  level  of  scattered  field  was  much  larger  compared  to  the  desired  one.  The 
realized  metamaterial  structures  were  quite  complicated,  and  the  lack  of  “really  good” 
experimental  results  has  raised  the  question  where  are  the  problems  -  in  the  non-perfect 
general  electromagnetic  solver,  in  the  non-perfect  realization  of  a  very  complex  structure,  or 
in  some  hidden  (physical)  limitations  that  are  not  so  obvious.  Therefore,  there  is  a  need  to 
understand  the  real  essence  of  the  problem. 
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The  purpose  of  this  report  is  to  give  a  contribution  to  the  design  of  metamaterial  structures 
that  will  have  good  scattering  properties.  In  doing  so,  we  will  avoid  to  simply  press  the  button 
“optimization”  in  some  electromagnetic  solver.  We  will  try  to  develop  a  procedure  that  will 
start  with  the  physical  picture  of  how  our  structure  should  look  like  (in  the  cloak  case  first 
design  will  be  made  through  transformation  optics),  and  in  several  steps  we  will  arrive  to  the 
actual  realization.  It  should  be  noted  that  it  is  extremely  difficult  to  design  curved  periodic 
structures,  i.e.  it  is  much  easier  to  design  planar  periodic  structure  (in  this  case,  it  is  only 
necessary  to  design  a  unit  cell).  Therefore,  the  proposed  procedure  will  first  characterize  an 
equivalent  planar  structure,  its  realization  using  planar  metamaterials,  and  in  the  last  step  will 
include  the  effect  of  curving  the  metamaterials  needed  to  obtain  the  desired  geometry. 
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2  PROJECT  OBJECTIVE 
AND  REALIZED  OUTCOMES 
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Project  objective  and  realized  outcomes 


The  work  in  the  project  has  been  divided  in  the  following  tasks: 


•  Development  of  the  method  -  systematic  approach  for  characterizing  metamaterial 
structures.  The  method  should  use  results  of  test  cases  calculated  by  a  general 
electromagnetic  solver,  and  the  selected  test  cases  will  depend  on  the  geometry  of  the 
considered  structure.  Two  types  of  parameters  will  result  from  the  analysis:  (a) 
reflection  and  transmission  coefficients  of  the  considered  structure,  and  (b) 
permeability  and  permittivity  material  tensors  of  the  considered  periodic  structure 
(homogenization  approach  -  needed  for  designing  the  first  prototype  where  the 
desired  constitutive  parameters  are  determined  by  the  program  that  can  analyze 
multilayer  anisotropic  homogeneous  structure).  Furthermore,  the  rigorous  full-wave 
parameters  of  the  considered  structure  will  be  also  determined. 

•  Estimation  of  parameters  of  planar  multilayer  periodic  structures.  The  considered 
structures  will  be  periodic  in  two  dimensions,  and  it  will  have  finite  number  of  layers 
(all  the  layers  will  be  the  same).  Two  basic  configurations  will  be  investigated  - 
configurations  in  free-space  and  configurations  that  are  part  of  some  transmission  line 
(e.g.  a  planar  or  cylindrical  waveguide).  It  will  be  discussed  how  the  number  of  layers 
influences  the  constitutive  parameters  of  the  structure.  Furthermore,  it  will  be 
investigated  is  it  possible  to  treat  each  layer  separately,  or  is  it  necessary  to  analyze  the 
whole  structure  (which  is  much  more  complex  electromagnetic  problem). 

•  Estimation  of  parameters  of  curved  multilayer  periodic  structures.  It  will  be 
investigated  how  the  radius  of  curvature  influences  the  constitutive  parameters.  In 
particular,  the  cylindrical  waves  will  be  analyzed  in  details. 


The  main  realized  outcomes  of  the  project  are: 

•  We  have  proposed  the  procedure  for  designing  metamaterial  structures.  The  procedure 
consists  of  four  steps.  In  the  first  step  a  device  is  designed  using  some  physical 
method  (like  transformation  optics)  and  then  it  is  realized  using  layers  with 
homogeneous  anisotropic  materials.  In  the  second  step  an  equivalent  planar  multilayer 
structure  is  characterized  using  reflection  &  transmission  coefficient  method.  In  the 
third  step  a  metamaterial  periodic  structure  is  designed  (using  some  general 
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electromagnetic  solver)  which  has  the  same  scattering  properties  like  multilayer 
homogeneous  structure  from  the  second  step.  In  the  last  step  the  designed 
metamaterial  structure  is  curved,  and  the  curvature  effects  are  taken  into  account  using 
e.g.  conformal  mapping  approach. 

•  We  have  investigated  properties  of  different  procedures  for  extracting  effective 
constitutive  parameters.  It  was  shown  that  the  method  is  quite  accurate  (and  the 
extracted  effective  parameters  quite  useful  in  the  design  procedure)  if  the  metamaterial 
elements  are  not  resonant.  However,  for  frequencies  around  resonant  the  structure 
becomes  narrow-banded,  and  it  is  almost  impossible  to  predict  its  behavior. 

•  We  have  made  a  detailed  investigation  of  the  procedure  that  connects  the  considered 
curved  structure  and  the  equivalent  planar  structure.  It  is  shown  that  with  the  method 
based  on  conformal  mapping  it  is  possible  to  accurately  predict  the  curvature  effects. 
Furthermore,  with  the  Bloch  theory  for  cylindrical  structures  it  is  possible  to 
efficiently  analyze  multilayer  cylindrical  periodic  structures. 

•  We  have  given  several  design  examples  of  metamaterial  structures  containing  different 
types  of  periodic-cell  elements.  In  all  cases  the  results  indicate  that  the  proposed 
design  procedure  is  quite  accurate  and  enables  fast  design  of  the  desired  metamaterial 
structure. 

•  We  have  developed  several  algorithms  and  programs  needed  for  characterization  of 
metamaterial  structures.  First,  we  have  developed  an  extension  of  the  G1DMULT 
algorithm,  which  calculates  Green’s  functions  of  planar,  cylindrical  and  spherical 
multilayer  structures,  to  calculate  Green’s  functions  of  anisotropic  multilayer 
structures  as  well.  Theoretical  formulation  is  given  in  the  Appendix,  and  the  algorithm 
was  used  in  designing  “ideal”  metamaterial  structures  built  from  homogeneous 
anisotropic  layers.  Furthermore,  we  have  developed  moment-method  based  programs 
for  analyzing  planar  and  cylindrical  periodic  structures.  Their  description  is  also  given 
in  the  Appendix. 


We  sincerely  hope  that  this  report  will  help  researchers  in  understanding  where  are  the 
critical  points  when  designing  metamaterial  structures,  and  will  result  in  an  efficient  and 
successful  designing  procedure. 
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3  DESIGN  OF  METAMATERIAL  STRUCTURES 
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Design  of  metamaterial  structures 


3.1.  EXAMPLE:  CLOAK  DESIGN  USING  TRANSFORMATION  OPTICS 

The  possibility  of  cloaking  objects  using  a  metamaterial  cover  has  extensively  been  studied 
in  the  last  several  years  (see  e.g.  [9]-[15]).  In  the  metamaterial  cloak  approach,  material  has 
been  used  to  render  a  volume  effectively  invisible  to  incident  radiation,  i.e.  to  squeeze  space 
from  a  volume  into  a  shell  surrounding  the  concealment  volume.  Coordinate  transformations 
that  are  used  for  cloak  design  do  not  influence  the  form  of  Maxwell’s  equations,  but  they 
affect  permittivity  and  permeability  tensors  (s  and  p,  respectively),  making  the  needed 
materials  spatially  varying  and  anisotropic.  When  viewed  externally,  the  concealed  volume 
and  the  cloak  both  appear  to  have  the  propagation  properties  of  free  space,  i.e.  they  appear 
invisible  to  electromagnetic  waves. 


Figure  3.1.  A  sketch  of  the  analyzed  structure. 


As  an  example  we  consider  a  circular-cylindrical  cloak  oriented  with  the  axis  in  z- 
direction  as  shown  in  Fig.  3.1.  The  inner  and  outer  radii  of  the  cloak  are  denoted  by  a  and  b. 
The  permeability  and  permittivity  tensors  in  the  cylindrical  coordinate  system  are  given  with 
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Figure  3.2.  (a)  a  field  line  in  free  space  with  the  Cartesian  coordinate  grid  in  the  background, 
(b)  the  distorted  field  line  with  the  background  coordinates  distorted  in  the  same  manner. 


For  the  cloak  design,  a  coordinate  transformation  which  compresses  free  space  from  the 
cylindrical  region  0  <  r  <  b  into  the  concentric  cylindrical  shell  a  <  r  ’  <  b  is  applied,  where  a 
and  b  represent  the  cloak  inner  and  outer  radius,  respectively  (the  constitutive  structure 
parameters  are  considered  in  the  cylindrical  coordinate  system).  Presuming  that 
electromagnetically  inhomogeneous  metamaterials  offer  more  freedom  in  design  than  regular 
materials  found  in  nature,  it  is  possible  to  imagine  for  electromagnetic  fields  to  be  controlled 
and  redirected  at  will.  The  design  of  a  cloaking  structure  is  based  on  the  transformation  of 
space,  i.e.  the  stretching  and  the  pulling  of  a  Cartesian  coordinate  system  as  shown  in  Figure 
3.2.  The  Maxwell’s  equations  are  invariant  under  spatial  transformations,  i.e.  they  take  the 
same  form  in  the  transformed  coordinate  system: 


V'xE'=  -jcoju'H',  V'xH'=  jcoe'  E' 


(3.3) 


The  permeability  and  permittivity  tensors  (1  and  £  are  related  to  the  tensors  ft  and  £ 
in  the  original  space  via  equations: 


P  = 


ApA7 


£ =  ■ 


AeA1 


det  A  det  A 

where  A  is  the  Jacobi  matrix  of  the  spatial  transformation 

dx'/dx  dx'/dy  dx'/dz 


A  = 


dy'/dx  dy'/dy  dy'/dz 
dz '/  dx  dz '/  dy  dz '/  8z 


(3.4) 


(3.5) 
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In  the  considered  case  the  circular  annular  cloak  is  obtained  using  the  following 
transformation: 


,  b-a 

p=a  +  — — p,  (p  =  <p,  z  =  z, 

b 

and  the  Jacobi  matrix  is  given  by 
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(3.6) 


(3.7) 


This  transformation  allow  us  to  design  a  2D  cylindrical  electromagnetic  invisibility  cloak 
by  compressing  all  the  fields  from  the  region  p  <  b  into  the  region  a  <  p  <  b  so  that  the 
permittivity  and  permeability  of  the  medium  located  within  the  cloaking  shell  are  free  to  take 
any  value  without  contributing  to  electromagnetic  scattering.  For  a  <  p  <  b  the  permittivity 
and  permeability  values  are  [9]: 


Spp  Ppp 


p-a 

P 


£  — 


P 


p-a 


p. 


b-a 


p-a 

P 


(3.8) 


Note  that  all  components  of  the  tensors  are  functions  of  radius,  which  implies  a  very 
complicated  metamaterial  design.  Such  a  cloak  is  referred  to  as  an  ideal  cloak  that  has  not  yet 
been  realized.  However,  several  cloak  designs  have  been  reported  which  claim  to  work 
properly  when  illuminated  with  a  normal  incident  plane  wave  of  specific  polarization  [16], 
[17],  that  is,  with  either  the  electric  field  or  the  magnetic  field  parallel  to  the  z-axis  (TMZ 
polarization  or  TEZ  polarization,  respectively). 


The  metamaterial  cloak,  designed  by  Schurig  et  al.  [16],  is  intended  for  use  with  TMZ 
polarization  of  the  incident  wave  at  microwave  frequencies.  It  is  clear  that  only  szz,  pf)f,  and  p^ 
components  are  relevant  in  this  case.  Therefore,  it  is  possible  to  vary  the  values  of 
permeability  and  permittivity  tensors  as  long  as  the  products  szz-ppp  and  szz-p^  are  kept  the 
same.  A  significantly  simplified  design  can  be  achieved  by  fixing  values  of  szz  and  p^  (e--  is 
determined  by  supporting  dielectric  material  and  p^  is  simply  equal  to  1),  and  letting  only  ppp 
to  vary  along  the  radial  direction: 


Pm  = 


fp-<P 
P 


=  1 


(3.9) 


£zz 


b-a 
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The  metamaterial  reported  in  [16]  used  the  unit  cell  of  approximately  cubic  shape,  filled  with 
one  split-ring  resonator.  The  realized  cloak  consisted  of  ten  concentric  cylinders,  each  of 
which  was  three  unit  cells  tall  (see  Figure  3.3). 


(a)  (b) 

Figure  3.3.  (a)  sketch  of  the  first  experimentally  developed  cloak,  (b)calculated  and  measured 

scattered  field  (Schurig  et  al,  2006  [16]). 


The  described  cloak  was  prototyped  and  embedded  into  scattering  chamber.  The 
distribution  of  the  vertical  component  of  the  electric  field  around  the  cloak  was  scanned,  and 
the  obtained  results  showed  some  decrease  of  the  scattered  field  (see  Fig.  3.3.b).  Although 
this  experiment  proved  the  basic  idea  of  cloaking,  the  level  of  scattered  field  has  not  been 
quantified. 

Similarly  to  the  realization  of  a  so-called  TMZ  cloak,  there  is  another  realization  of 
metamaterial  cloak  by  Cai  et  al.  [17].  This  cloak  is  intended  to  work  with  TEZ  polarized  wave 
at  optical  frequencies.  In  TEZ  cloak  only  components  pzz,  epp  and  are  relevant,  thus  it  is 
again  possible  to  simplify  the  metamaterial  design  by  allowing  only  epp  to  vary  in  radial 
direction: 


£ pp 


Saa.  = 


f  u  p-a V 


\b-a  j 
f  u 


P  J 


\b  —  ay 

Pzz=l- 


(3.10) 


Again,  the  products  epp'pzz  and  are  kept  the  same.  The  required  distribution  of  epp  can 

be  realized  using  metal  wires  embedded  in  a  dielectric  material,  i.e.  with  a  thin  wire  plasma¬ 
like  metamaterial.  Here,  the  effective  permittivity  varies  in  radial  direction  and  its  real  part 
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has  to  exhibit  the  required  behavior  (0  <  Refe^}  <  1)  with  negligible  imaginary  part  (i.e  with 
negligible  losses)  [17].  First  practical  realization  in  microwave  frequency  range  has  used 
split-rings  to  obtain  required  permittivity  distribution  (see  Fig.  3.4  and  references  [18]  and 
[19]  for  details). 


Figure  3.4.  Sketch  of  the  experimentally  developed  TEZ  cloak  ([18],  [19]) 


(a) 


(c) 
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Figure  3.5.  Experimental  demonstration  of  the  TEz  cloak. 
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Table  3.1  gives  the  comparison  of  the  permeability  and  permittivity  values  of  the  ideal, 
TMZ  and  TEZ  cloaks.  It  should  be  noted  that  for  some  of  the  proposed  cloak  designs  there  is  no 
need  to  use  metamaterial  structures.  For  example,  hard  surfaces  have  been  used  in  the  design 
of  supporting  struts  of  reflector  antenna  feeds  [6],  and  parallel-plate  waveguide  structures 
have  been  used  to  reduce  scattering  from  dielectric  or  metallic  obstacles  [21],  [22], 


Table  3.1.  The  values  of  permeability  and  permittivity  tensors  of  the  considered  cloaks 

(a  AND  b  REPRESENT  THE  INNER  AND  OUTER  RADIUS  OF  THE  CLOAK.  AND  P  IS  THE  RADIAL 

COORDINATE) 


Ideal  cloak  [9] 

TM  cloak  [16] 

TE  cloak  [17] 

p-a 

PfP  = 

f  ^ 

p-a 

2 

„  J  b  Yf p~a Y 

^  pp  ^pp 
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v  P  J 
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_  Hdxb  ~ 
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£  -f  M2 

\b-a) 

f  b  Y  p-a 

r»  —  //  —  ' 

szz 

{  b  ^ 

2 

Pzz  =  1 

^  ZZ  r'zz  7 

V b-aj  p 

yb-a  j 

In  last  several  years  different  versions  of  cylindrical  cloaks  are  experimentally  realized, 
both  in  microwave  and  optical  frequency  ranges.  Unfortunately,  for  some  realizations  (like 
Schurig  cloak  [16]  or  TEZ  cloak  [19])  there  are  no  published  data  in  open  scientific  literature 
which  will  give  information  about  the  level  of  scattered  field  or  about  the  obtained  frequency 
bandwidth.  In  other  words,  these  experiments  have  only  proved  the  basic  idea  of  cloaking. 

There  is  only  one  experimentally  realized  metamaterial-based  cloak  for  which  data  about 
the  obtained  invisibility  and  bandwidth  are  available  [20].  It  was  realized  at  Duke  University 
and  the  cloak  is  nearly  identical  to  the  cloak  realized  by  Schurig  et  al  [16].  The  central 
frequency  was  10  GHz,  and  the  total  scattering  width  (of  the  bare  PEC  cylinder)  was  reduced 
in  the  frequency  range  9.91  to  10.14  GHz.  In  other  words,  the  obtained  bandwidth  was  230 
MHz  or  2.3%  (which  corresponds  well  to  theoretical  predictions).  The  reduction  of  total 
scattering  width  was  24%  (i.e.  the  invisibility  gain  was  only  1.32),  and  from  the  measurement 
report  [20]  it  is  not  clear  if  such  small  obtained  invisibility  gain  was  consequence  of  losses  or 
imperfections  in  the  experimental  realization.  The  experimentally  measured  total  scattering 
width  is  given  in  Figure  3.6. 
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Figure  3.6.  The  measured  total  scattering  width  cT  of  an  experimentally  realized  TMZ  cloak 
[20].  For  comparison,  aj  of  an  uncloaked  metal  cylinder  is  also  given. 

The  experimental  results  give  rise  to  a  natural  question,  why  the  realized  cloak  has  such 
poor  performance.  There  are  two  basic  reasons  for  that  (additionally  to  the  problem  of 
tolerances  and  of  the  accuracy  of  electromagnetic  solvers),  as  described  below.  In  previous 
work  we  have  developed  a  program  that  analyses  multilayer  uniaxial  cylinders  (see 
[EOARD]),  which  enabled  us  to  analyze  the  “ideal”  realization  of  the  Schurig  cloak.  We  have 
selected  the  same  dimensions  and  the  working  frequency  as  in  the  experimental  model 
realized  by  Schurig  et  al.  [16]  (see  Tables  3.1  and  3.2). 

Table  3.2.  Dimensions  of  the  cylindrical  cloak 


dimension 

(cm) 

Dimension 

(wavelength) 

Inner  diameter  2 a 

5.42 

1.44 

Outer  diameter  2b 

1 1.78 

3.34 

Thickness  of  one  layer 
(10-layer  case) 

0.318 

0.09 

First  we  have  investigated  the  influence  of  the  stepwise  approximation  of  the  continuous 
variation  of  constitutive  parameters.  The  required  has  been  approximated  by  1-  to  10-  steps 
piecewise  constant  functions  in  order  to  analyze  how  the  subtlety  of  the  approximation  of 
radial  anisotropy  influences  the  total  scattering  width,  i.e.  “the  invisibility”.  The  normalized 
total  scattering  width  is  given  in  Fig.  3.7,  and  it  can  be  seen  that  for  structures  with  more  than 
5  layers  there  is  practically  no  improvement  of  the  obtained  invisibility  -  the  obtained  total 
scattering  width  reduction  is  around  3.  The  main  reason  for  such  a  small  gain  was  the 
reflection  of  the  incident  wave  from  the  cloak  surface  due  to  the  impedance  mismatch.  The 
level  of  scattered  power  can  be  approximated  by 
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Here  cloak  and  pec  denotes  the  total  scattering  width  of  the  cloaked  and  uncloaked 
cylinder,  respectively.  If  we  compare  the  calculated  value  (0.15)  and  approximated  value 
(0.09)  we  can  conclude  that  the  impedance  mismatch  is  the  main  reason  why  Schurig  cloak 
always  scattered  large  amount  of  power.  The  angular  variation  of  scattered  electric-field  and 
the  field  distribution  in  the  vicinity  of  the  cloaked  object  is  given  in  Fig.  3.8.  and  it  can  be 
seen  that  the  reflection  from  the  cloak  surface  causes  ripples  in  the  electric-field  distribution. 


Figure  3.7.  Normalized  total  scattering  width  vs.  number  of  layers  for  TMZ  cloak 

(TMZ  polarization). 


(a)  (b) 

Figure  3.8.  (a)  Normalized  bistatic  scattering  width  (TMZ  and  TEZ  polarizations)  and  (b) 
electric  field  distribution  in  the  vicinity  of  the  cloaked  object  (TMZ  polarization)  of  10-layers 
realization  of  TMZ  cloak. 
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Next  we  have  investigated  the  bandwidth  of  the  cloak.  The  frequency  dependence  of  the 
magnetic  permeability  (realized  with  some  kind  of  split  rings)  is  given  by  the  so-called 
Lorentz  model  [23]: 


M#=1- 


f  2  -  f 2 

J  mp  J  0 
/2-/o2-i>/' 


(3.12) 


Here /is  the  frequency  of  the  signal,/,,,,  denotes  the  frequency  at  which  jueJf=0  for  the  lossless 
case  (null-point  of  the  function),  fo  is  the  frequency  at  which  fieff  diverges  (the  pole  of  the 
function),  while  the  factor  y  represents  the  losses. 

In  an  attempt  to  approximate  the  experimental  Schurig  cloak  we  have  to  select  one  free 
parameter.  The  bandwidth  results  of  the  TMZ  cloak  strongly  depend  on  the  ratio  between  the 
frequencies  fmp  (the  frequency  at  which  jueJf=  0  for  the  lossless  case)  and  fo  (the  frequency  at 
which  neff  diverges  for  the  lossless  case).  Since  this  ratio  depends  on  the  particular  realization 
of  the  split-ring  resonator,  we  have  investigated  the  dependence  of  the  cloak  bandwidth  on 
that  ratio.  In  Figure  3.9  we  have  illustrated  this  dependence.  It  can  be  seen  that  the  bandwidth 
is  between  0.24%  and  2.4  %,  depending  on  this  ratio.  Therefore,  one  can  conclude  that  if  the 
cloak  is  made  from  resonant  type  of  elements  the  bandwidth  will  be  quite  narrow,  typically 
several  percent  (as  it  was  demonstrated  in  [20]).  However,  if  we  decide  to  select  non-resonant 
elements  (like  wires  that  are  used  in  TEZ  cloak  realization),  the  bandwidth  can  be  much  larger 
(at  least  one  order  of  magnitude). 


Figure  3.9.  Normalized  total  scattering  width  as  a  function  of  ratio  fmp  I  fo  for  TMZ  cloak 

(TMZ  polarization). 
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3.2.  EXTRACTION  OF  PARAMETERS 

Previous  work  in  invisibility  gives  rise  to  several  questions: 

1.  Why  the  developed  cloak  with  the  reduced  variation  of  constitutive  parameters  (i.e. 
cloaks  made  from  uniaxial  materials)  scatters  the  electromagnetic  field,  even  in  the 
ideal  case? 

2.  Is  it  possible  to  realize  a  cloak  with  the  reduced  variation  of  constitutive  parameters 
that  is  entirely  invisible  (at  least  at  the  central  frequency)? 

3.  Is  it  possible  to  realize  a  cloak  with  the  reduced  variation  of  constitutive  parameters 
that  works  for  oblique  direction  of  incident  wave? 

4.  What  is  the  bandwidth  of  a  cloak  that  is  built  using  metamaterial  layers?  Is  it  possible 
to  enlarge  the  bandwidth  of  such  a  cloak? 

5.  Is  it  possible  to  realize  a  cloak  that  is  invisible  to  pulse  excitation  (i.e.  to  radar)? 

We  believe  (and  we  will  try  to  demonstrate  that  in  this  report)  that  with  smart  designing 
procedure  it  is  possible  to  avoid  most  of  possible  failures,  i.e.  designing  decisions  that  will 
result  in  a  developed  structure  with  poor  electromagnetic  performance  (the  explanation  of  the 
phenomena  resulting  in  the  first  question  is  already  given  in  section  3.1). 

Before  starting  we  should  note  that  metamaterial  structures  are  extremely  hard  to  design 
since  they  have  a  lot  of  small  details  (small  compared  to  the  wavelength)  which  results  in  a 
large  computation  time  when  analyzing  such  structures  using  general  electromagnetic  solvers. 
Therefore,  it  is  not  possible  simply  to  connect  the  electromagnetic  solver  with  a  optimization 
routine  and  press  the  button  “Optimize”. 

All  these  lead  to  conclusion  that  the  first  design  should  be  made  using  (anisotropic) 
homogeneous  materials.  This  could  be  done  using  e.g.  transformation  optics  (as  described  in 
section  3.1)  or  simply  by  connecting  the  fast  electromagnetic  solver  with  the  optimization 
routine  (homogeneous  structures  can  be  analyzed  very  quickly).  We  believe  (due  to 
complexity  of  possible  realizations)  that  the  design  obtained  in  such  a  way  is  actually  the  best 
possible  design.  In  other  words,  the  practical  realization  will  always  have  worse  performance 
compared  to  this  ideal  model. 

In  order  to  connect  the  ideal  design  (based  on  homogeneous  multilayer  anisotropic 
structure)  and  “real”  metamaterial  structure  we  should  somehow  connect  the  considered 
periodic  structure  with  effective  constitutive  parameters,  i.e.  with  parameters  that  describe  the 
scattering  properties  of  the  considered  structure. 

Extraction  of  parameters  or  homogenization  is  a  procedure  in  which  complex  periodic 
structure  is  represented  with  a  simpler  structure.  Usually  that  simpler  structure  can  be 
described  with  several  numbers,  i.e.  constitutive  parameters.  Procedure  deals  with  equations 
in  heterogeneous  materials,  with  a  periodic  pattern,  in  which  dimensions  of  a  unit  cell  tend  to 
zero.  Method  is  based  on  the  consideration  of  two  length  scales  associated  with  the 
microscopic  and  macroscopic  phenomena.  Constitutive  parameters  of  homogeneous  materials 
do  not  depend  on  spatial  variable.  On  the  other  hand,  if  the  material  is  not  homogeneous, 
parameters  of  the  material  effectively  depend  on  the  position  in  space.  Parameters  of  material 
with  periodic  pattern  are  periodic  functions  of  space  variables.  In  certain  cases  the  period  is 
very  small  with  respect  to  other  dimensions  appearing  in  the  problem.  In  those  cases  solution 
is  approximately  the  same  as  the  solution  for  a  homogeneous  material  with  corresponding 
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constitutive  parameters  [24].  The  homogenization  procedure  in  electromagnetic  theory 
introduces  constitutive  parameters  of  a  discrete  media  as  a  result  of  averaging  of  the  Maxwell 
equations  for  microscopic  fields.  The  procedure  makes  sense  if  its  result  can  be  used  for 
solving  boundary  problems  of  electrodynamics  where  a  discrete  set  of  scattering  particles  are 
replaced  by  a  sample  of  continuous  media  [25].  The  necessity  for  a  good  homogenization 
procedure  is  related  with  the  request  for  a  fast  electromagnetic  solver,  i.e.  within  procedure  of 
metamaterial  design. 

The  constitutive  parameters  can  be  obtained  by  numerical  calculations  or  by  some 
approximate  analytical  model.  In  numerically  based  methods  we  simply  calculate  the  ratio  of 
electromagnetic  field  components.  Approximate  analytical  methods  are  based  on  geometrical 
properties  of  periodic  structure  so  there  are  some  practical  problems  with  implementing  this 
for  complicated  metamaterial  structures.  Furthermore,  due  to  complex  geometrical  pattern, 
some  unit  cells  of  periodic  structures  exhibit  bianisotropic  behavior. 


3.2.1  Extraction  based  on  reflection  and  transmission  coefficient 

In  paper  [26]  the  method  for  extraction  of  constitutive  parameters  from  measured  S 
parameters  is  presented.  The  idea  is  to  extract  constitutive  material  parameters  from  reflection 
and  transmission  coefficients  for  a  unit  cell  in  all  three  spatial  directions.  In  all  those 
directions  our  unit  cell  is  illuminated  by  a  normal  incident  plane  wave.  We  have  six 
incidences  on  a  unit  cell  (see  Figure  3.10).  Three  of  them  are  TE,  and  three  are  TM  modes.  S 
parameters  are  defined  in  terms  of  the  electric  and  magnetic  fields  for  the  TE  and  TM 
incidences.  In  every  direction  unit  cell  can  display  bianisotropic  properties  or  it  behaves  as 
isotropic  material.  In  directions  for  which  the  unit  cell  behaves  as  isotropic  material  we  can 
use  formulas  from  [27],  For  all  these  directions  the  retrieval  equations  are  the  same. 


Figure  3.10.  S-matrix  approach  for  determining  constitutive  parameters.  The  method  is  based 
on  measuring  or  calculating  Sn  and  S21  for  six  test  cases. 
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Elements  of  the  S  matrix  can  be  found  from  the  elements  of  the  T  matrix  as  follows.  The 
transfer  matrix  for  a  section  of  a  transmission  line  is  defined  as: 


F  =  T  F 


(3.13a) 


jco/uH 


(3.13b) 


cos  (nk0d) - sin  (nk0d) 

k0 

kn 

—  sin  (nk0d)  cos  (nk0d) 


(3.13c) 


In  [28]  there  is  a  connection  between  S  and  T  matrices: 


Tu  -  T22  +  (jk0T12  -  -f-) 
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(3.14.c) 


(3.14.d) 


For  a  homogeneous  slab  Tn=T22  and  det(7)=l.  Furthermore,  S  matrix  is  symmetric: 


-(—-jknTu) 

2  0  121 

Tn  +  2«>ri2+JL: 


(3.15) 
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Using  the  analytic  expression  for  T  matrix  elements  we  get: 
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^11  — 


S2i  = 


S22  =  ^(77  ~  Z)sin(nkQd) 
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(3.16) 


j  X 

cos  (nkQd)  —  -(—  +  Z)sm(nkQd ) 


Here  ko  denotes  the  wave  number  of  the  incident  wave  in  free  space,  d  is  the  thickness  of  the 
slab,  n  is  the  refractive  index,  R  is  the  half  space  reflection  coefficient.  From  these  two 
equations  we  can  get  expressions  for  index  of  refraction  n  and  characteristic  impedance  Z: 
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21 


Note  that  signs  in  equations  (3.17)  are  determined  from  conditions  for  passive  media. 
Effective  parameters  are  calculated  from: 


jj,  =  nZ 


(3. 18. a) 


(3.18.b) 


Table  3.3  describes  the  proposed  method  with  the  six  test  cases.  It  can  be  seen  that  each 
component  of  the  ft  and  £  tensors  is  calculated  twice.  Comparison  of  two  cases  gives  an 
indication  how  accurate  the  constitutive  parameters  are  calculated  (in  particular,  this  is  critical 
around  resonant  frequencies  of  metamaterial  elements). 

If  in  some  direction  material  behaves  as  a  bianisotropic  material  then  other  methods  have 
to  be  used.  Some  of  those  methods  are  explained  in  [29].  Those  formulas  are  not  enough  to 
determine  constitutive  parameters  because  in  some  situations  we  will  get  two  roots  of  some 
parameters  and  from  theory  we  cannot  conclude  which  one  is  the  physical  one. 
Recommendation  from  the  paper  is  to  use  one  which  gives  better  match  between  calculated 
and  simulated  or  measured  S  parameters. 


Table  3.3.  Extraction  of  parameters  by  S-matrix  approach  (six  test  cases  included). 


Case 

Z 

n 

Solved  components 

TE1 

4»*ley 

TM1 

My>£X 

TE2 

^y/£z 

^y£z 

Vy,£z 

TM2 

^-J£y 

^-£y 

Vz’£y 

TE3 

hlz£. 

sz,£x 

TM3 

MX’£Z 
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3.2.2  Extraction  based  on  Bloch  theory 

In  this  it  is  assumed  that  the  interactions  among  adjacent  frequency  selective  surfaces  (in 
the  normal  direction)  are  well  described  by  the  dominant  Floquet  mode,  and  that  other  modes 
have  negligible  coupling  effects.  From  this  assumption  equivalent  transmission  line  model  is 
adopted  in  which  every  FSS  layer  is  represented  by  an  equivalent  admittance  matrix.  In 
general,  the  admittance  of  FSS  layer  depends  on  frequency  and  on  direction  of  arrival  of 
electromagnetic  wave.  In  order  to  determine  the  FSS  admittance  the  full-wave  analysis  is 
adopted.  Isolated  FSS  is  illuminated  by  a  plane  wave  with  a  given  frequency  and  incident 
direction.  If  propagation  vector  is  in  the  principle  plane  of  the  FSS,  the  equivalent  admittance 
is  a  diagonal  matrix.  That  means  that  interaction  of  the  plane  wave  and  FSS  could  be  modeled 
with  decoupled  equivalent  transmission  lines,  one  for  TM  and  one  for  TE  polarization.  The 
whole  device,  which  is  made  as  a  multilayer  structure,  is  modeled  as  a  couple  of  periodically 
loaded  transmission  lines  that  are  analyzed  with  Bloch  theory  [30],  From  that  procedure  we 
get  propagation  constant  and  characteristic  impedance  of  the  supported  modes.  Then,  a 
homogeneous  medium  supporting  the  same  scattered  field  distribution  is  defined. 

One  of  the  problems  with  the  described  procedure  is  that  the  solution  is  not  unique.  Since 
the  transmission  line  model  only  provides  information  on  the  transverse  field  components, 
there  are  infinite  media  that  satisfy  the  boundary  conditions  of  the  solution.  All  those  media 
provide  the  same  dispersion  relation  and  scattering  coefficients  of  the  fundamental  modes  of 
the  periodic  pattern  but  they  are  characterized  by  different  constitutive  parameters,  so  all 
those  solutions  are  externally  equivalent.  In  order  to  uniquely  determine  the  effective 
constitutive  parameters  it  is  necessary  to  find  additional  information  on  the  longitudinal  field 
components  inside  periodic  structure.  It  is  assumed  that  electromagnetic  field  inside  material 
is  very  well  approximated  by  the  dominant  Bloch  mode  only.  Under  this  assumption,  a  closed 
form  expression  of  the  longitudinal  field  can  be  obtained  as  a  function  of  the  FSS  equivalent 
admittances.  That  step  provides  the  complete  characterization  of  the  homogenized  material. 
The  largest  weakness  of  this  method  lays  in  the  assumption  that  only  dominant  mode  is 
needed.  That  is  not  true  for  FSS  layers  which  are  close  to  each  other  because  in  that  situation 
we  need  more  modes,  i.e.  we  need  to  consider  evanescent  modes  as  well. 

In  the  first  step  we  calculate  the  reflection  coefficient  of  a  considered  FSS  using  some 
general  electromagnetic  solver.  Equivalently,  we  consider  the  network  from  Fig.3.11. 


Figure  3.11.  Equivalent  network  of  FSS  in  free-space. 

From  reflection  coefficient  T  we  calculate  the  normalized  FSS  admittance: 
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VTE,TM 

2fss 
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(3.18) 


The  Bloch  theory,  described  in  [30],  gives  the  parameters  of  the  equivalent  transmission  line: 


yTE^TM 

cos  (k™’™d)  =  cos  (kZod)  +  — sin  (kZQd) 


(3. 19. a) 


rvTE.TM 

L B  _  ?TE,TM 

rvTETM  ~  ^ B 
Z0 


2 sin (k^d)  -  jY^™  cos(kZ(d)  +  jV/Jf™ 
2 sin (kZQd)  -  jY^™  cos (kZQd)  -  jYpss™ 


(3.19.b) 


From  procedure  described  in  [31]  we  get: 
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(3.21) 


where  kp  is  a  propagation  constant  in  p  direction.  As  mentioned  above,  z  components  of 

constitutive  parameters  are  not  determined  in  this  procedure.  To  determine  them  we  need 
additional  information,  longitudinal  components  of  the  electromagnetic  fields.  The 
longitudinal  fields  are  constructed  by  properly  averaging  the  corresponding  components  of 
the  microscopic  fields  [32],  The  final  expressions  for  z-components  of  constitutive  parameters 
are: 
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(3. 23. a) 


(3.23.b) 


where  h  and  e  are  microscopic  fields.  With  these  two  expressions  we  have  obtained 
expressions  for  all  components  of  the  constitutive  parameters.  The  problems  with  practical 
implementation  lie  mostly  in  evaluation  of  microscopic  fields,  i.e.  the  accuracy  of  the  Bloch 
method  mostly  depends  on  the  accuracy  of  calculating  microscopic  near-fields. 
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3.3.  PROCEDURE  FOR  DESIGNING  METAMATERIAL  STRUCTURES 

When  designing  novel  type  of  structures  that  contain  metamaterial  layers  it  is  quite  tedious 
and  complicate  to  arrive  from  basic  idea  (obtained  for  example  using  transformation  optics)  to 
„real“  structure  that  will  be  produced  in  some  workshop.  The  main  obstacle  is  that  these  novel 
structures  (in  the  basic  design)  contain  layers  of  materials  that  have  permeability  or 
permittivity  smaller  than  one,  i.e.  such  materials  that  are  not  possible  to  find  in  nature. 
Therefore,  we  need  to  “cheat“  the  nature,  i.e.  to  use  periodic  structures  with  period  much 
smaller  than  the  wavelength,  i.e.  we  need  to  implement  the  layers  of  metamaterials. 
Therefore,  the  straightforward  path  is  to  design  a  structure  which  will  have  effective 
constitutive  parameters  equal  to  the  designed  ones.  However,  this  is  simple  only  for  some 
structures.  At  least  two  problems  appear  when  one  starts  to  design  structures  using  this 
approach.  In  more  details,  in  the  first  design  (i.e.  in  the  design  that  reflects  the  ideal  structure 
obtained  e.g.  with  ray  approach)  one  would  like  to  obtain  material  with  specific  either 
permeability  or  permittivity.  However,  the  practical  realization  will  also  introduce  variation  of 
other  quantity,  the  one  that  is  usually  kept  constant.  The  larger  problems  appear  from  the 
resonance  nature  of  elements  used  in  metamaterial  realization.  It  is  shown  [29]  that  around 
resonances  there  are  frequency  intervals  in  which  it  is  not  possible  to  accurately  define 
effective  constitutive  parameters,  i.e.  parameters  that  will  accurately  predict  the  behavior  of 
the  electromagnetic  waves  that  travel  through  such  a  structure. 

So,  the  natural  question  arise  that  if  is  not  possible  to  accurately  characterize  metamaterials 
using  effective  constitutive  parameters,  how  metamaterial  structures  should  be  designed.  Most 
often,  the  ideal  structure  is  described  using  (anisotropic)  permeability  and  permittivity,  so  we 
should  not  completely  forget  tensors  £  and  p.  But  what  we  actually  consider  is  how  the  waves 
are  propagating  through  such  a  structure.  Therefore,  the  description  that  is  more  practical 
from  the  design  point  of  view  is  description  of  transmitted  and  reflected  electromagnetic 
waves. 

We  will  describe  this  principle  on  one  example,  i.e.  how  to  design  circular-cylindrical 
cloak  (Figure  3.13).  In  the  ideal  design  cylindrical  layers  with  specific  permeability  and 
permittivity  are  obtained  using  transformation  optics  (see  previous  section).  These  layers  we 
would  like  to  obtain  using  cylindrical  periodic  structures.  From  CPU  time  and  memory 
requirements  it  is  much  faster  and  easier  to  analyze  planar  periodic  structures.  Therefore,  we 
will  actually  design  planar  multilayer  metamaterial  structure  and  then  we  will  band  it.  In 
designing  planar  structure  we  will  start  with  ideal  homogeneous  multilayer  structure  with 
desired  (anisotropic)  permeability  and  permittivity  and  then  we  will  calculate  reflection  and 
transmission  coefficients  for  set  of  incidence  angles  of  incidence  (and  for  polarization  of 
incoming  wave  of  interest).  This  dependency  of  reflection  and  transmission  coefficient  will  be 
our  cost  function  in  optimization  procedure,  i.e.  we  will  search  for  a  periodic  structure  that 
have  the  same  response  to  incoming  waves  as  the  desired  ideal  structure.  This  is  now  not  a 
difficult  task  since  we  need  to  analyze  only  one  periodic  cell  (multilayer  in  perpendicular 
direction),  which  is  quite  fast  for  all  nowadays  commercial  electromagnetic  solvers.  This 
designing  procedure  is  illustrated  in  Figure  3.12. 
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Design  of  the  structure  using  anisotropic 
homogeneous  materials  approach 
(ideal  design) 
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Electromagnetic  characterization  of 
equivalent  planar  homogeneous  structure 
(planar  version  of  idelal  structure) 


Design  of  planar  metamaterial  structure 
that  has  similar  electromagnetic  properties 
as  planar  homogenouus  structure 


Modification  of  metamaterial  structure 
—*■  new  design  has  included 
curvature  effects 


Figure  3.12.  The  proposed  designing  procedure  of  metamaterial  structures 
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Figure  3.13.  Illustration  of  the  proposed  designing  procedure:  circular-cylindrical  cloak 
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3.3.1  Characterization  of  a  planar  metamaterial  layer 

As  a  first  step  let  us  characterize  planar  multilayer  anisotropic  structures.  The  “real“ 
structures  usually  are  not  planar,  but  this  first  step  in  designing  procedure  will  give  accurate 
starting  point  for  final  design,  i.e.  for  final  optimization  with  general  electromagnetic  solver. 
For  example,  the  cylindrical  cloak  is  actually  circular-cylindrical  multilayer  anisotropic 
structure.  The  second  step  in  the  design  process  (Fig.  3.12)  solves  the  sub-problem  how  to 
realize  equivalent  planar  multilayer  structure  using  metamaterial  approach. 

From  mathematical  point  of  view  the  question  is  how  to  characterize  planar  multilayer 
structure  and  to  use  this  characterization  as  a  cost  function  in  optimizing  metamaterial 
structure.  The  original  structure  is  shown  in  Figure  3.14.  Since  lots  of  applications  are 
actually  scattering  electromagnetic  problem,  we  should  also  characterize  the  considered 
structure  using  scattering  approach.  For  example,  if  we  would  like  to  characterize  the 
structure  in  Fig.  3.14  that  contains  PEC  ground  plane,  there  are  several  possibilities  of 
defining  the  cost  function  /cost.  Some  of  them  are  (the  reflection  coefficient  of  a  homogeneous 
anisotropic  layer  and  of  the  periodic  metamaterial  structure  is  denoted  by  /lom  and  rma, 
respectively): 

1.  Structure  with  PEC  ground  plane  — »  phase  of  reflection  coefficient  as  a  function  of 
incoming  angle.  The  cost  function  is  equal: 

U  =  |>g(r„((2))-arg(rloJ0,))|  (3,24) 

i'=l 

2.  Structure  without  PEC  ground  plane  — ►  amplitude  of  reflection  coefficient  as  a 
function  of  incoming  angle.  The  cost  function  is  equal: 

f„„  =  S|r„,  W,  )|  -  |rl0„(6(  )|  (3.25) 

(=1 

3.  Structure  without  PEC  ground  plane  —*■  amplitude  and  phase  of  reflection  and 
transmission  coefficients  as  a  function  of  incoming  angle. 

4«,  =Elr'™,W)-rh„J^)|  (3.26) 

;= l 


In  order  to  illustrate  the  accuracy  of  these  methods  let  us  consider  the  following  examples. 
The  first  one  is  simple  one-layer  structure,  of  thickness  2  cm  and  with  ground  plane,  which 
has  permeability  equal  to  1.0,  and  uniaxial  permittivity  sx  =  ez  =  1.0  and  sy  =  0.5  (see  Figure 
3.14).  The  working  frequency  is  4  GHz.  Such  a  structure  can  be  realized  using  metal  strips, 
for  example  with  periodicity  Px  =  3  cm.  In  order  to  get  an  estimate  of  the  strip  width,  we  have 
calculated  the  effective  constitutive  parameters  of  a  dielectric  slab  (sr  =  1 .0  in  the  considered 
case)  containing  the  periodic  strips  in  the  middle  of  the  structure  (see  Fig.  3.15).  For  that 
purpose  we  have  developed  a  Moment  Method  code  for  planar  periodic  structure  (we  have 
named  this  code  pfss,  and  its  description  is  given  in  the  Appendix  of  this  report).  From 
calculated  reflection  and  transmission  coefficient  for  normal  incidence  we  have  calculated  the 
constitutive  parameters  using  the  method  explained  in  section  3.2.  The  calculated  results  show 
that  for  obtaining  sy  =  0.5  we  need  to  select  the  strip  width  W  =  0.9  mm.  The  results  also  show 
that  the  strips  introduce  change  in  effective  permeability,  i.e.  for  the  considered  case  the  x- 
component  of  the  p  tensor  is  pix  =  1.145. 
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Figure  3.14.  Dielectric  slab  with  periodic  metal  strips 


Figure  3.15.  Effective  constitutive  parameters  of  a  dielectric  slab  with  periodic  metal  strips 


Figure  3.16.  Phase  of  reflection  coefficient  of  a  grounded  anisotropic  dielectric  layer  and  of 
metamaterial  realization.  Anisotropic  layer  is  characterized  with  sx  =  sz  =  1.0,  sy  =  0.5  and 
px  =  1.0.  The  used  cost  function  in  optimization  routine  is  defined  by  eq.  (3.24).  Parameters  of 
strips  are  listed  in  Table  3.4. 
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Figure  3.17.  Phase  of  reflection  coefficient  of  a  grounded  anisotropic  dielectric  layer  and  of 
metamaterial  realization.  Anisotropic  layer  is  characterized  with  sx  =  sz  =  1.0,  sy  =  0.5  and 
px=  1.146.  The  used  cost  function  in  optimization  routine  is  defined  by  eq.  (3.24).  Parameters 
of  strips  are  listed  in  Table  3.4. 


Figure  3.18.  Amplitude  and  phase  of  reflection  coefficient  of  an  anisotropic  dielectric  layer 
and  of  metamaterial  realization.  Anisotropic  layer  is  characterized  with  sx  =  sz  =  1.0,  sy  = 
0.5  and  with  px  =  1.0.  The  used  cost  function  in  optimization  routine  is  defined  by  eq.  (3.25). 
Parameters  of  strips  are  listed  in  Table  3.4. 
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Figure  3.19.  Amplitude  and  phase  of  reflection  coefficient  of  an  anisotropic  dielectric  layer 
and  of  metamaterial  realization.  Anisotropic  layer  is  characterized  with  sx  =  e2  =  1 .0,  sy  = 
0.5  and  with  px  =  1.146.  The  used  cost  function  in  optimization  routine  is  defined  by  eq. 
(3.25).  Parameters  of  strips  are  listed  in  Table  3.4. 


Figure  3.20.  Amplitude  and  phase  of  reflection  coefficient  of  an  anisotropic  dielectric  layer 
and  of  metamaterial  realization.  Anisotropic  layer  is  characterized  with  sx  =  e2  =  1 .0,  sy  = 
0.5  and  with  px  =  1.0.  The  used  cost  function  in  optimization  routine  is  defined  by  eq.  (3.26). 
Parameters  of  strips  are  listed  in  Table  3.4. 
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Figure  3.21.  Amplitude  and  phase  of  reflection  coefficient  of  an  anisotropic  dielectric  layer 
and  of  metamaterial  realization.  Anisotropic  layer  is  characterized  with  sx  =  ez  =  1 .0,  sy  = 
0.5  and  with  px  =  1.146.  The  used  cost  function  in  optimization  routine  is  defined  by  eq. 
(3.26).  Parameters  of  strips  are  listed  in  Table  3.4. 


Figure  3.22.  Amplitude  and  phase  of  reflection  coefficient  of  an  anisotropic  dielectric  layer 
and  of  metamaterial  realization.  Anisotropic  layer  is  characterized  with  sx  =  ez  =  1 .0,  sy  = 
0.5  and  with  px  =  1.146.  The  width  of  the  strips  is  determined  using  effective  constitutive 
parameters  approach. 
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Table  3.4.  Comparison  of  different  realizations  of  strip-grid  layer  obtained  using 

VARIOUS  OPTIMIZATION  APPROACHES  (CONSIDERED  COST  FUNCTIONS  ARE  DEFINED  WITH 

EOS.  (3.12) -(3.14)). 


Method 

Cost 

function 

Px  of 

homogeneous 

anisotropic 

structure 

Optimum  width 
of  the  strips 
(mm) 

Value  of 
effective  ey 
(Figure  3.14) 

Value  of 
cost  function 

1 

(3.12) 

1.0 

3.002 

0.292 

0.35560 

1 

(3.12) 

1.146 

0.2 

0.534 

0.15631 

2 

(3.13) 

1.0 

1.356 

0.444 

0.83416 

2 

(3.13) 

1.146 

0.986 

0.485 

0.09344 

3 

(3.14) 

1.0 

0.979 

0.487 

0.47885 

3 

(3.14) 

1.146 

1.092 

0.474 

0.07079 

effective 

- 

1.0 

0.9 

0.5 

0.89252 

effective 

- 

1.146 

0.9 

0.5 

0.11154 

Figures  3.16  -  3.21  show  the  comparison  of  amplitude  and  phase  of  reflection  coefficient 
of  metamaterial  structure  obtained  using  three  different  approaches  (in  the  first  case  with  a 
PEC  ground  plane,  so  the  amplitude  of  reflection  coefficient  is  equal  to  one).  Two  types  of 
quality  control  can  be  examined:  similarity  between  reflection  coefficients  (as  a  function  of 
incidence  angle)  and  effective  constitutive  parameters  of  the  considered  metamaterial 
structure  (when  the  considered  metamaterial  structure  is  not  resonant  the  effective  parameter 
approach  gives  accurate  results).  For  that  reason,  for  each  considered  case  we  have  listed  the 
optimized  parameters,  the  value  of  the  cost  function  and  the  value  of  effective  sy,  see  Table 
3.4.  It  can  be  seen  that  the  first  approach  (phase  of  the  reflection  coefficient  from  the 
grounded  metamaterial/anisotropic  layer)  gives  quite  good  results  at  the  first  glance  (see 
Figures  3.16  and  3.17).  However,  the  optimized  width  of  the  strips  is  quite  different  in 
comparison  to  the  value  obtained  by  the  effective  parameters  approach.  Therefore,  the  first 
approach  does  not  give  accurate  results  (since  the  presence  of  the  PEC  plane  “kills“ 
information  about  the  structure),  so  it  is  much  better  to  take  out  the  PEC  ground  plane 
(although  it  is  present)  and  to  consider  the  reflection  coefficient  of  the  structure  without  the 
ground  plane. 

The  second  approach  (no  PEC  ground  plane  present  and  only  the  amplitude  of  reflection 
coefficient  is  considered  in  the  definition  of  cost  function)  gives  better  results,  i.e.  more 
similar  to  the  effective  parameters  approach.  However,  in  this  case  the  influence  of  the 
induced  permeability  (px  ^  1.0)  is  clearly  seen.  If  we  state  that  px  =  1.0  (Fig.  3.18)  the 
agreement  between  the  amplitude  of  the  reflected  signal  from  two  structures  is  not  good.  One 
can  drastically  improve  the  accuracy  of  the  metamaterial  realization  simply  by  adding  small 
permeability  to  the  structure  (Fig.  3.19).  However,  in  “real”  practical  situations  it  is  not 
obvious  if  we  are  allowed  to  simply  add  this  small  permeability  to  the  designed  structure. 

The  last  approach  is  if  the  used  cost  function  takes  into  account  both  the  amplitude  and  the 
phase  of  the  structure  (Figures  3.20  and  3.21  and  Table  3.4).  It  can  be  seen  that  the  accuracy 
of  this  method  is  the  largest  (although  not  so  different  comparing  to  the  previous  case)  and 
that  the  optimized  width  of  the  strips  is  the  closest  to  value  we  would  expect.  Therefore,  in  the 
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rest  of  the  report  we  will  use  the  cost  function  defined  by  eq.  (3.26)  in  the  optimization 
process. 

Final  comparison  gives  the  results  for  the  case  when  we  would  use  just  effective 
constitutive  parameters  approach,  see  Figure  3.22  and  Table  3.4.  Two  values  of  effective  px 
are  considered:  px  =  1.0  and  px  =  1.146.  The  width  of  the  strips  is  determined  using  effective 
constitutive  parameters  approach,  i.e.  the  width  W  =  0.9  mm  is  taken  from  the  Figure  3.15.  It 
can  be  seen  that  there  is  good  matching  for  normal  incidence  (almost  perfect  matching  if  we 
take  into  account  that  effective  permeability  is  different  from  free-space  for  strip-loaded 
substrates),  and  that  this  matching  is  deteriorated  for  angles  different  from  6  =  0.  If  we 
calculate  the  error  of  this  approach  (i.e.  the  cost  function  using  in  the  optimization  procedure), 
it  can  be  seen  that  this  value  is  larger  comparing  to  the  results  obtained  by  optimizing  the 
strips  width  value,  see  Table  3.4. 


3.4.  CURVATURE  EFFECTS 


In  the  previous  section  we  have  described  a  design  procedure  in  which  the  main  step  is 
metamaterial  realization  of  the  equivalent  planar  multilayer  anisotropic  structure.  The  real 
device  will  contain  curved  structures,  so  we  need  to  transform  back  the  planar  structure.  Two 
possibilities  will  be  discussed  in  this  section:  conformal  mapping  and  Bloch  theory  for 
cylindrical  structures. 

3.4.1.  Conformal  mapping 

In  order  to  take  curvature  effects  into  account  we  could  use  techniques  that  actually  we 
have  used  for  designing  metamaterial  devices  -  transformation  optics.  However,  this 
procedure  can  lead  to  complicated  mathematical  expressions  that  we  could  evaluate  only 
numerically.  Therefore,  let  us  consider  a  subgroup  of  transformation  techniques  that  are 
suitable  for  this  task  -  conformal  mapping. 

A  function  defined  in  the  complex  plane  f(z )  =  f(x  +  jy )  =  u(x,  y)  +  jv(x,  y)  is  conformal 
if  it  satisfies  the  requirements  for  an  analytic  function,  i.e.  Caushy-Riemann  conditions: 


du  dv  ,  dv  du 

—  =  —  and  —  = - . 

dx  dy  dx  dy 


(3.27) 


On  the  other  side,  transformation  of  the  area  in  the  free-space  is  given  by  equation  (3.4): 


p  =  £  = 


AAr 
det  A 


(3.28) 


The  two-dimensional  conformal  mapping  can  be  expressed  in  three-dimensions  as  [33]: 
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(u,v,z)  =  (u(x,y),v(x,y),z)  (3.29) 

By  taking  into  account  the  Caushy-Riemann  conditions  (3.27)  one  can  calculate  Jacobian 
matrix: 


du/dx 

du/dy 

du/dz 

du/dx 

du/dy 

0" 

A  = 

dv/dx 

dv/dy 

dv/dz 

= 

-du/dy 

du/dx 

0 

dz/dx 

dz/dy 

dz/dz 

0 

0 

1 

The  determinant  IAI  is  equal: 

|  A|  =  ( du/dxf  +  (du/dy)2  =  ( du/dxf  +  ( dv/dxf 
Therefore,  the  constitutive  parameters  of  the  transformed  medium  are  equal: 


p  =  £  = 


AAr 
det  A 


1  0  0 

0  1  0 

0  0  1/1 A I 


(3.30) 


(3.31) 


(3.32) 


In  order  to  understand  properties  of  the  obtained  transformation  let  us  consider  an  example 
-  transformation  of  planar  layer  with  strips  into  a  circular-cylindrical  layer  (see  Figure  3.23). 
Let  us  suppose  that  we  have  already  designed  the  metamaterial  realization  of  anisotropic  layer 
with  permittivity  smaller  than  one.  Now  we  would  like  to  see  how  this  structure  is 
transformed  to  a  cylindrical  layer. 


Figure  3.23.  Transformation  of  planar  to  cylindrical  metamaterial  structure. 


The  conformal  mapping  that  supports  this  transformation  is  given  by 
c Z=x  +  jy\  x 

P strips  y  P strip sfi)  * 

f(z)  =  exp(z/ps7r.  J  =  exp([A(n/Jn(p)  +  jpsln,j]/  Pslrips)=  p •  (cos^+ ,/sin </>) .  (3.33) 

The  factor  /Strips  is  introduced  to  “match”  transformation  from  the  arc  length  (i.e.  width  of  the 
strips)  to  the  angular  coordinate.  In  other  words,  the  real  part  of  z  complex  variable  represents 
the  thickness  of  the  considered  layer,  and  the  imaginary  part  represents  the  angular 
coordinate.  The  determinant  of  matrix  A  is  equal: 
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|A|(p,  (j>)  -  (i du/dxf  +  (du/dy)2 


—j—  (exp(x/ ps,rips)f  •  (cos2  </>  +  sm2</>)=  -f— .  (3.34) 

Pstrips  P strips 


For  example,  let  us  consider  the  same  strip-grid  structure  like  in  the  previous  section 
(Px  =  3  mm,  W  =  0.9  mm,  d  =  2  cm,  /  =  4  GHz).  This  structure  correspond  to  an  anistropic 
layer  with  sx  =  sz=  1.0,  sy  =  0.5  and  px  =  1.146.  We  would  like  to  consider  three  cylindrical 
structure  radii  -  around  0.5  Ao,  1 .0  Ao  and  1.5  Ao  (i.e.  the  number  of  strips  in  the  azimuthal 
direction  were  8,  16  and  24,  respectively).  Figure  3.24.  shows  these  three  cases,  i.e.  we 
compared  the  scattered  field  from  a  homogeneous  anisotropic  cylindrical  layer  and  from  the 
strip  structure  obtained  by  the  conformal  mapping.  Note  that  only  the  cylindrical  structure 
geometry  is  obtained  using  conformal  mapping,  and  that  the  scattering  properties  of  periodic 
cylindrical  structure  is  calculated  by  in-house  developed  code  for  the  analysis  of  cylindrical 
periodic  structures  (cylindrical  FSS;  description  of  the  code  is  given  in  the  Appendix).  The 
results  show  excellent  agreement  between  scattered  field  for  structures  with  larger  radius,  and 
some  discrepancy  in  the  back-scattered  field  for  the  structures  with  small  radius. 


Figure  3.24.  Scattered  field  from  the  periodic  strip-grid  structure  and  from  the  homogeneous 
anisotropic  cylindrical  shell,  (a)  pstnps  =  0.5 A0,  (b)  pstrips  =  1.0A0 ,  (c)  pstrips  =  l.5A0 . 
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3.4.2.  Bloch  theory  for  cylindrical  structures 

The  Bloch  theory  of  periodic  transmission  lines  can  be  extended  to  the  cylindrical  structures. 
The  basic  block  of  the  method  is  a  unit  cell  that  is  periodically  repeated  (see  Figure  3.24).  In 
our  case  we  have  quasi-periodical  structure  in  the  radial  direction  and  periodical  structure  in 
the  (f)  -direction. 


Figure  3.25.  Decomposition  of  periodic  unit  cell  into  3  parts. 

The  considered  unit  cell  will  be  decomposed  into  3  parts,  and  for  each  part  the  ABCD  matrix 
will  be  determined.  Then  the  ABCD  matrix  of  the  unit  cell  is  equal  to  the  product  (cascade)  of 
ABCD  matrices  of  each  part. 


A 


p=0 


Figure  3.26.  First  segment  of  the  unit  cell. 

The  ABCD  matrices  describe  the  considered  system  in  the  following  way: 


un 

I 


A  B 
C  D 


U 


n+1 


n+l 


(3.35) 


Here  subscript  n  represents  position  on  the  radial  transmission  line,  i.e.  n  denotes  the  radial 


position  pn .  Situation  on  the  line  is  described  with  linear  combination  of  two  independent 

solutions  of  the  Helmholtz  equation.  In  our  situation  these  solutions  describes  incoming  and 
outgoing  wave,  so  representation  matrix  for  a  line  is: 
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(3.36) 


Here  An  and  Bn  represents  the  E-  and  H-field,  and  Zo  is  the  characteristic  impedance  of  the 
free-space  (or  of  the  homogeneous  media  that  fills  the  transmission  line).  From  those  two 
equations  we  can  deduce  equation  for  voltage  and  current  on  the  line: 


\u  1 

\u  ^  1 

n 

=  NM-1 

n+ 1 

I 

I 

n 

n+ 1 

(3.37) 


which  is  equivalent  to  the  ABCD  representation,  so  the  ABCD  matrix  is  equal  to  NM -1 .  For 
inverse  computation  we  need  to  calculate  the  determinant  of  the  matrix  M  .  For  easier 
computation  we  introduce: 


det  M  = 


1 


X  =  kpPn 


y  =  kpPn+ 1 


1  -4  i 

JZ0  7rkpPn+ 1 


(3.38) 

(3.39) 


Using  this  information  we  get: 


M~x 


P^pPn+l 

4 


-jZoHfXy) 

jZ0H^(y) 


(3.40) 


Now  we  have  all  needed  expressions 
the  cylindrical  coordinate  system. 


A  B 
C  D 


=  NM-1  = 


m- 


X 


H(,!\ 


X 


to  determine  the  ABCD  matrix  of  transmission  line  in 


j^y 

'  to)' 

JZ0  K  ’  \ 

4 

jiry 

4 


-  H^(y)'H^\x) 


-jZ»  <SS\x)H{;\y)  -  H'P(y)H^(x)  > 


4-  <gi'\x)'H^(y)'  -  H^(y)'H^(x)'  )  -  C'M'-fff  to)  -  H(»(y)H«\x)' 
jZq 


(3.41) 
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We  can  verify  this  matrix  via  the  analogy  with  planar  case.  For  large  argument  of  the  Hankel 
functions  (i.e.  when  cylindrical  structure  is  almost  planar)  we  get: 


A  B 
C  D 


pry 

4 


—=  r'V.  j)  -  — —  e*e~*j  jZ0 

■Kjxy  n-sfxy 

— - —  jejxe  jy(  j) - —  jejye  jx(  j)  — . 

iZ0  n^fxy  n^fxy  irjxy 


n^fxy 


-  eJxe~JV 


H)  - 


ir^fxy 

2 


TT^fxy 


ejye~jxj 


-i 


cos{x  -  y )  jZ0  sin(x  -  y) 

j 

—  sin  (x  —  y)  cos  (x  —  y) 

Z0 


(3.42) 


So  for  large  radii  the  ABCD  matrix  of  the  cylindrical  transmission  is  equivalent  to  the  ABCD 
matrix  of  the  transmission  line  in  the  Cartesian  coordinate  system,  i.e.  we  have  only  phase 

shift  present  (and  small  adjustment  of  the  amplitude  since  cylindrical  waves  decay  as  l/  Jp  ). 


If  the  FSS  admittance  is  known  it  is  easy  to  compute  the  ABCD  matrix.  It  is  equal  to  the 
admittance  ABCD  matrix  in  the  Cartesian  coordinate  system: 


A 

B 

1 

0 

C 

D 

y 

1  FSS 

1 

(3.43) 


ABCD  matrix  of  the  unit  periodic  cell  is  equal  to  the  product  of  ABCD  matrices  of  unit  cell 
parts.  To  make  expressions  easier  to  code  we  will  introduce  slightly  different  notation: 

1  =  V» 

V  =  k  p  J  (3.44) 

n-\ — 

2 

z  =  k  p  ,, 

p>  n+ 1 

From  Figure  3.25,  derivations  from  previous  two  paragraphs  and  using  notations  (3.44)  we 
derive: 


u  ' 

n 
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Pr  n  +  0.5 
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-  H*\x)'H*\y)  -  H^(y)H[2\x)' 

1  i 

n  +  ~ 

L  2  -1 

u 

1 

n  +- 

2 

I 

1 

n  +- 
2 


1 

Y 

FSS 


0 

1 


u  / 
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(3.45) 


Distribution  A:  Approved  for  public  release;  distribution  is  unlimited. 


Constitutive  Parameters  of  Metamaterial  Structures  Used  for  Invisible  Cloak  Realization  42 


To  introduce  even  simpler  notation  let  us  introduce  the  following  abbreviations: 

V  =  H J 

H?  z'  =  Hi 

zrl  Tt2  TT  1  TT 2  _  TT 

x  y  y  x  xy 

zrl  r?2  TT  1  tt 2  _  TT 

x  y  y  x  xy 

fjl  fj2  zrl  fj2  _  TT 

x  y  y  x  xy 

The  ABCD  matrix  of  the  unit  cell  in  cylindrical  coordinate  system  is  equal: 


(3.46) 


A  B' 

■  \2 
J7V 

C  D 

,  4  , 

yz 


HxyHyz  ^O^FSS^xy^yz  ^xy^yz 

—  H  H  +YFqqH  H  +—H  H 

■  ry  xy  yz  FSS  yx  yz  ry  yx  yz 

^Zn  Zn 


~J^oHxyHyz  ZQYFSSHxyHyz  jZ0HxyHzy 
jZ,YFssHyxHyz  d-  HzyHyx 


TT  7T 

11  xy 11  yz 


(3.47) 


After  grouping  terms  that  represent  elements  of  the  ABCD  matrix: 

AL  =  Hjyz  ~  jZ0YFSSHxyHyz  -  HxyHyz 
BL  =  ~jZ0HxyHyz  -  ZlYFSSHxyHyz  -  jZQHxyHzy 

CL  =  +  YfssHvHv  +  pH  vHtz 

DL  =  II  H  -  jZ0YFSSHyxHyz  +  HzyHyx 

we  get  the  ABCD  matrix  in  simplified  notation: 


X' 

/  .  \ 

jir 

2 

AL 

BL 

yn+ 1 

4 

yz 

CL 

DL 

^n+ i 

Using  the  Bloch  Theorem  in  cylindrical  coordinate  system  we  can  conclude  that  ratio  of 
eigenvalues  behaves  as: 


lMp„+1)l=  [aT 

\Mpn)\  \pn+ 1 


(3.50) 


We  can  test  that  conclusion  on  one  example.  We  can  use  SRR  unit  cell  with  d= 8  mm  from  the 
paper  [29]  and  for  validation  we  have  calculated  eigenvalues  of  the  above  ABCD  matrix,  see 
Figure  3.27.  We  took  the  distance  1.6  cm  from  the  origin  as  arefemce  point  and  we 
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investigate  the  relative  change  of  the  eigenvalues  of  the  ABCD  matrix  (we  always  compare 
eigenvalues  for  two  neighboring  /^-coordinates).  If  we  go  further  away  from  the  origin,  ratio 
will  become  more  and  more  accurate  because  of  the  large-argument  approximation. 


Figure  3.27.  Validation  of  expression  (3.50).  Distance  1.6  cm  from  the  origin  is  taken  as  a 

reference  point. 


The  described  ABCD  matrix  approach  represent  a  powerful  tool  for  analyzing  multilayer 
periodic  cylindrical  structures.  One  just  need  to  analyze  each  layer  separately  (more  precisely 
one  just  needs  to  determine  YFSS  of  each  cylindrical  layer  with  FSS),  and  the  ABCD 
approach  construct  the  solution  for  the  multilayer  case.  This  approach  also  simplifies  the 
optimization  process  since  it  is  possible  to  separate  the  original  complex  problem  into  several 
simpler  sub-problems. 
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4  EXAMPLES 
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EXAMPLES 


4.1  PERMITTIVITY  TAILORING  -  LAYERS  WITH  STRIPS 

The  basic  structure  that  enables  permittivity  values  smaller  than  one  (or  even  negative)  are 
periodic  wires  or  strips  (their  shape  can  be  curved  in  order  to  obtain  the  desired 
performances).  We  would  like  to  test  the  described  design  procedure  on  a  more  complex 
structure  (comparing  the  one  given  in  the  previous  section).  Let  us  consider  a  four-layer 
homogeneous  anisotropic  structure  with  the  following  parameters  (the  working  frequency  is  4 
GHz): 

Table  4.1.  Constitutive  parameters  of  anisotropic  multilayer  homogeneous 

STRUCTURE. 


layer 

thickness  (cm) 

permittivity 

1 

2 

£x=£z=l-0,  sv  =  0.125 

2 

2 

£x  =  £z  =  1.0,  £v  =  0.  25 

3 

2 

£x=£z=1.0,  Sy  =  0.3755 

4 

2 

£x  =  £z  =  1.0,  £y  =  0.5 

We  would  like  to  realize  such  a  structure  using  layers  with  strip-grids  (see  Fig.  4.1).  Let  the 
periodicity  of  the  strips  is  Px  =  3  cm  (like  in  the  previous  example).  The  parameter  that  we 
would  like  to  determine  is  the  width  of  the  strips  in  each  layer  (therefore,  we  have  an 
optimization  problem  with  four  unknowns).  Two  approaches  are  possible:  using  the  design 
procedure  that  we  have  described  in  previous  section,  or  simply  to  take  the  width  that  is 
determined  for  each  layer  by  considering  the  effective  parameters  of  isolated  layer  (i.e.  the 
width  values  are  determined  from  Figure  3.15).  The  results  obtained  by  two  design  methods 
are  described  in  Figures  4.2  and  4.3.  In  the  optimization  procedure  there  are  two  possible 
approaches  regarding  the  permeability  values  of  the  equivalent  homogeneous  anisotropic 
structure:  either  to  simply  state  that  px  =  1.0  (Figure  4.2)  or  to  take  expected  value  of  the 
induced  x-component  of  the  permeability  tensor  px  =  1.2  (Figure  4.3).  It  can  be  seen  that  if  we 
include  in  the  optimization  procedure  the  expected  induced  value  of  the  permeability  tensor 
the  obtained  scattering  parameters  are  much  more  close  to  the  scattering  parameters  of  the 
homogeneous  anisotropic  multilayer  structure.  Furthermore,  the  determined  widths  of  the 
strips  are  much  more  close  to  the  values  obtained  by  considering  effective  constitutive 
parameters  of  a  single-layer  structure  (i.e.  the  width  is  determined  from  Fig.  3.15).  One 
should  also  notice  that  if  we  select  px  =  1.0  the  widths  of  the  strips  are  not  monotonically 
decreasing,  i.e.  the  obtained  structure  does  not  have  physical  picture  like  the  case  when  px  = 
1.2.  These  results  are  summarized  in  Table  4.2. 


Distribution  A:  Approved  for  public  release;  distribution  is  unlimited. 


Constitutive  Parameters  of  Metamaterial  Structures  Used  for  Invisible  Cloak  Realization  46 


W 


Figure  4.1.  Multilayer  metamaterial  structure  containing  periodic  metal  strips 


(a) 


Figure  4.2.  Amplitude  and  phase  of  reflection  coefficient  of  an  anisotropic  dielectric 
multilayer  structure  and  of  metamaterial  realization.  Anisotropic  layer  is  characterized  with 
constitutive  parameters  given  in  Table  4.1  (permeability  px  =  1.0).  The  width  of  the  strips  is 
given  in  Table  4.2. 
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(a)  (b) 

Figure  4.3.  Amplitude  and  phase  of  reflection  coefficient  of  an  anisotropic  dielectric 
multilayer  structure  and  of  metamaterial  realization.  Anisotropic  layer  is  characterized  with 
constitutive  parameters  given  in  Table  4.1  (permeability  px  =  1.2).  The  width  of  the  strips  is 
given  in  Table  4.2. 


Table  4.2.  Comparison  of  different  realizations  of  strip-grid  layer  obtained  using 

OPTIMIZATION  APPROACH  AND  EFFECTIVE  PARAMETERS  APPROACH. 


Method 

Px  of 

homogeneous 

anisotropic 

structure 

Optimum 
width  of  the 
strips  (mm) 

Value  of 
effective  sy 
(Figure  3.4) 

Value  of 
cost  function 

optimization 

method 

1.0 

Wi  =  2.79 

W2  =  5.76 

W3  =  2.01 

W4  =  2.41 

£y  =  0.310 
sy  =  0.075 
£y  =  0.379 
£y  =  0.343 

0.60819 

optimization 

method 

1.2 

Wi  =  6.36 

W2  =  3.00 

W3  =  2.84 

W4  =  0.67 

£y  =  0.029 
sy  =  0.292 
sy  =  0.306 

Sy  =  0.536 

0.23724 

effective 

parameters 

1.0 

Wi  =  5.11 

W2  =  3.53 

W3  =  2.05 

W4  =  0.9 

£y  =  0.125 

Sy  =  0.25 
sy  =  0.375 

Sy  =  0.5 

4.33642 

effective 

parameters 

1.2 

Wi  =  5.11 

W2  =  3.53 

W3  =  2.05 

W4  =  0.9 

£y  =  0.125 

Sy  =  0.25 
sy  =  0.375 

Sy  =  0.5 

0.72906 

Here  we  will  mention  just  one  possible  application  of  such  structures.  In  some  cases  we 
would  like  to  make  the  radiation  of  basic  antennas  (like  dipoles)  more  directive.  That  can  be 
made  by  placing  a  layer  that  has  permittivity  smaller  than  one  in  front  of  the  non-directive 
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antenna  (see  Fig.  4.4).  In  this  example  we  will  place  the  four-layer  metamaterial  structure 
from  Fig.  4.3  (and  Table  4.2)  in  front  of  a  dipole  antenna,  i.e.  this  multilayer  structure  will  act 
as  a  collimating  lens.  The  parameter  of  interest  is  radiation  pattern  in  H-plane,  see  Fig.  4.5.  It 
can  be  seen  that  we  managed  to  obtain  the  directive  pattern,  and  that  practically  there  is  no 
difference  in  results  obtained  using  homogeneous  anisotropic  structure  and  real  metamaterial 
periodic  structure.  Reflections  from  the  metamaterial  planar  lens  introduce  the  losses  of  2  dB 
in  the  main  beam  in  the  main  beam  direction,  but  this  can  be  improved  with  a  careful  design. 


Figure  4.4.  Dipole  antenna  with  planar  metamaterial  lens. 


5 


e  (degrees) 


Figure  4.5.  H-plane  radiation  pattern  of  a  dipole  antenna  with  planar  metamaterial  lens. 
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4.2  PERMEABILITY  TAILORING 

-  LAYERS  WITH  SPLIT-RING  RESONATORS 


Split-ring  resonators  (SRRs)  are  usually  used  when  permeability  with  negative  or  close  to 
zero  values  is  needed.  The  main  disadvantage  of  SRR  elements  is  their  resonant  behavior.  In 
other  words,  the  negative  and  near-to-zero  values  of  permeability  appear  for  frequency  in  the 
vicinity  of  the  resonant  frequency  (more  precisely,  for  frequencies  larger  than  resonant 
frequency).  This  resonant  behavior  makes  practical  designs  very  difficult.  As  it  is  shown  in 
previous  chapter  the  bandwidth  of  such  structures  is  very  narrow.  Even  if  we  are  „satisfied“ 
with  narrow  bandwidth  it  is  still  extremely  hard  to  design  some  practical  structure.  The  main 
reason  again  lies  in  the  resonant  behavior  and  non-perfectness  of  all  electromagnetic  solvers 
(including  the  commercial  ones).  This  non-perfectness  is  shown,  for  example,  in  the  shift  of 
calculated  resonant  frequency  (disadvantage  of  all  commercial  solvers),  which  is  very  critical 
when  designing  metamaterial  structures.  Furthermore,  due  to  non  symmetry  of  the  structure, 
bianisotropy  effects  are  induced,  i.e.  the  constitutive  relations  are  now 


D  =  £  E  +  ^-H 
B=^ E+p H 


This  resonant  and  bianisotropy  behavior  of  SRRs  makes  extraction  of  effective  parameters 
extremely  difficult.  With  the  described  method  of  extraction  from  the  reflection  and 
transmission  coefficients  of  normal  incident  plane  wave,  we  actually  extract  each  the 
parameter  twice  (see  Table  3.3).  The  agreement  between  the  extracted  parameters  for  two 
different  excitation  direction  indicate  us  what  is  the  accuracy  of  the  proposed  method.  For 
example,  let  us  try  to  extract  the  constitutive  parameters  from  the  SRR  structure  from  Figure 
4.6.  The  parameters  of  the  structure  are:  a  =  5.25  mm,  g  =  d  =  0.125  mm,  w  =  3.75  mm.  The 
most  important  components  of  the  constitutive  parameters  are  pyy  (the  negative  or  near-to- 
zero  permeabilty  that  we  want  to  obtain)  and  szz.  Note  that  Hy  is  inducing  an  electric  dipole  in 
the  z-direction,  and  Ez  is  inducing  a  magnetic  dipole  in  the  y-direction  (therefore,  these  two 
components  are  strongly  coupled).  Figures  4.7  -  4.9  shows  the  extracted  effective  parameters 
obtained  by  CST  Microwave  studio,  and  Figure  4.10  the  effective  permeability  obtained  by 
in-house  moment  method  program.  The  attempt  to  extract  the  effective  parameters  had  lead  to 
the  following  conclusions: 

(1)  The  method  of  extraction  of  parameters  is  not  accurate  close  to  the  resonant 
frequency..  It  can  be  seen  that  all  three  considered  components  of  effective  constitutive 
parametsres  (pyy,  szz.  and  sxx)  strongly  depends  on  the  method  of  extraction  around  the 
resonant  frequency  (i.e.  in  the  frequency  range  of  interest).  This  effect  is  mostly  due  to 
bianisotropy  nature  of  SRRs.  However,  even  if  we  take  bianisotropy  into  account  (see  [29]) 
the  retrieval  method  is  still  not  accurate  around  the  resonant  frequency. 

(2)  The  calculated  resonant  frequency  depends  on  the  selected  electromagnetic  solver.  In 

order  to  illustrate  that  we  have  analyzed  the  same  structure  with  the  in-house  developed 
electromagnetic  solver  based  on  the  Moment  Method.  Figure  4.10  shows  the  extracted  values 
of  pyy.  It  can  be  seen  that  the  extracted  curve  has  the  same  shape  like  for  the  results  obtained 
by  CST  Microwave  Studio,  but  the  resonant  frequency  is  more  than  1  GHz  lower.  This 
(unknown)  shift  of  the  resonant  frequency  is  disadvantage  of  all  nowadays  electromagnetic 
solvers,  and  only  the  measurements  can  show  how  large  is  that  difference. 
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(3)  There  are  strong  mutual  coupling  effects  between  the  layers  with  SRR  elements.  To 

illustrate  this  effect  let  us  consider  4-laer  structure  consisting  of  identical  SRR  elements  (the 
unit  cell  therefore  contains  4  cubes  from  the  Figure  4.11).  It  can  be  seen  that  the  mutual 
coupling  effects  are  extremely  strong  around  the  resonant  frequency,  and  that  it  is  not  possible 
to  say  what  are  the  values  of  the  equivalent  effective  parameters  for  frequencies  around  the 
resonant  frequency.  Such  strong  effects  are  not  present  when  building  layers  from  non¬ 
resonant  elements  (or  from  resonant  elements  with  different  enough  resonant  frequency).  To 
illustrate  that  we  have  considered  4-layer  structure  of  strip  grids  (of  identical  period  and  strip 
width;  the  period  Px  was  30  mm,  the  layer  thickness  d  was  30  mm  and  the  width  of  the  strips 
W  was  0.9  mm).  It  can  be  see  that  there  is  almost  no  difference  of  extracted  effective 
parameters  for  one-layer  and  four-layer  structures,  except  around  resonant  frequencies  of  a 
Fabry-Perot  resonator  built  from  partly  reflective  mirrors  (i.e.  partly  reflective  periodic 
structures).  Therefore,  in  practical  metamaterial  designs  the  resonant  elements  should  be 
avoided  (if  possible). 


Figure  4.6.  Sketch  of  the  unit  cell  of  split-ring  structure. 
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Figure  4.7.  Extracted  (xyy  for  the  structure  in  Fig.  4.6;  (a)  TE2  way,  (b)  TM1  way  of  retrieval 
of  effective  parameters  (definitions  are  given  in  Fig.  3.10  and  Table  3.3). 
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Figure  4.8.  Extracted  £xxfor  the  structure  in  Fig.  4.6;  (a)  TM1  way,  (b)  TE3  way  of  retrieval 
of  effective  parameters  (definitions  are  given  in  Fig.  3.10  and  Table  3.3). 


f  (GHz)  f  (GHz) 

(a)  (b) 


Figure  4.9.  Extracted  szz  for  the  structure  in  Fig.  4.6;  (a)  TE2  way,  (b)  TM3  way  of  retrieval 
of  effective  parameters  (definitions  are  given  in  Fig.  3.10  and  Table  3.3). 
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Figure  4.10.  Extracted  pyy  of  the  structure  in  Fig.  4.6  obtained  by  the  Moment  Method. 


Figure  4.11.  Extracted  pyy  of  the4-layer  structure  in  Fig.  4.6;  (a)  TE2  way,  (b)  TM1  way  of 
retrieval  of  effective  parameters  (definitions  are  given  in  Fig.  3.10  and  Table  3.3). 


Distribution  A:  Approved  for  public  release;  distribution  is  unlimited. 


Constitutive  Parameters  of  Metamaterial  Structures  Used  for  Invisible  Cloak  Realization  53 


Figure  4.12.  Extracted  pyy  of  the  strip-grid  4-layer  structure  in  Fig.  4.6;  (a)  TM1  way,  (b)  TE3 
way  of  retrieval  of  effective  parameters  (definitions  are  given  in  Fig.  3.10  and  Table  3.3). 
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4.3  ROTATORS  OF  POLARIZATION 

In  this  section  we  will  give  an  example  of  designing  metamaterial  device  that  cannot  be 
realized  using  effective  constitutive  parameters,  but  that  can  be  designed  using  reflection  and 
transmission  coefficients  approach. 

Rotators  of  polarization  and  polarizers  find  their  applications  in  microwave  systems,  where 
it  is  required  to  change  the  polarization  of  an  electromagnetic  wave.  It  is  very  well  known  that 
a  periodic  structure  of  parallel  wires  or  strips  can  rotate  plane  of  polarization  of  a  linearly 
polarized  wave  by  a  desired  angle  over  a  broad  frequency  range.  In  order  to  obtain  a  high 
accuracy  with  this  design  in  a  broad  frequency  range,  and  improve  the  transmission 
coefficient  at  the  same  time,  it  is  necessary  to  use  a  great  number  of  screens.  That 
automatically  means  significantly  increase  in  both  the  polarizer  thickness  and  weight  of  the 
system.  Utilizing  a  smaller  amount  of  screens  the  change  from  linear  to  elliptical  or  circular 
polarization  is  observed  in  higher  frequency  bands  [35]. 

Our  realizations  will  be  made  from  multilayered  thin  screens  that  are  two  dimensional 
periodic  structures  with  period  much  smaller  than  the  wavelength  (therefore  is  called 
metamaterial  structure).  Different  geometries  for  unit  cell  realizations  will  be  considered.  Unit 
cell  will  be  made  from  two  separated  surfaces.  One  will  be  at  the  edge,  and  second  will  be  in 
the  middle  of  the  cell.  Idea  is  to  discretely  rotate  the  middle  part  of  unit  cell  in  separated 
screens.  After  several  steps  screens  elements  will  be  rotated  for  90  degrees.  So,  at  the  front 
side  element  is  oriented  in  one  direction  and  on  the  back  side  element  is  oriented 
perpendicular  to  the  first  direction.  It  is  reasonable  to  expect  that  in  some  frequency  range  the 
designed  device  works  as  a  rotator.  All  realizations  are  simulated  in  CST  Microwave  Studio. 

1st  design 

As  a  first  design  we  considered  a  slot  in  the  middle  of  the  patch  which  is  surrounded  by  the 
free  space.  There  is  a  cascade  of  7  thin  screens  and  slot  is  rotated  by  15  degrees  between 
screens  in  negative  direction.  System  is  excited  by  a  plane  wave  orthogonally  to  the  surface  of 
the  system. 


Figure  4.13.  1st  design:  slot  in  a  patch,  a=1 0  mm,  L=8  mm,  R=4.5  mm,  b=lmm 
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This  unit  cell  doesn’t  have  a  load  (metal  insert)  in  the  slot  so  we  cannot  expect  resonant 
behavior  in  the  frequency  range  of  interest.  In  all  other  realizations  an  inverted  geometry  will 
be  used. 

2nd  design 

This  realization  is  a  complementary  structure  to  the  previous  design.  From  the  Babinet 
principle  we  can  expect  broader  bandwidth  (frequency  range  with  Sn  below  -10  dB).  Starting 
geometry  is  equal  to  the  starting  geometry  in  previous  realization  (a=10  mm,  L=8  mm,  R=4.5 
mm,  b=lmm). 


Figure  4.15.  2nd  design:  dipole  in  a  slot  surrounded  by  metal 
(a=  10  mm,  L=8  mm,  R=4.5  mm,  b=l  mm) 


SZmaxTE(0,0)ZmaxTE(0,0)dB 
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SZminTE(0,0)ZmaxTE(0,0)dB 


Figure  4.17.  S21  parameter  of  the  2nd  design,  coupling  between  TE  input  and  output  waves 


SZminTM(0,0)ZmaxTE(0,0)dB 


Figure  4.18.  S21  parameter  of  the  2nd  design,  coupling  between  TE  input  and  TM  output  wave 


First  geometry  have  bandwidth  of  approximately  1  GHz  (between  11.5  and  12.5  GHz),  and 
in  that  range  cross-coupling  S21  coefficient  (TE  input  plane  wave  and  TM  output  plane  wave) 
is  below  -2dB.  In  that  range  2nd  design  behaves  as  an  excellent  rotator. 

If  we  want  rotator  in  lower  frequency  range  we  have  to  increase  the  unit  cell  size.  Also,  we 
can  change  the  cell  geometry  (and  keep  the  same  outer  cell  dimensions),  and  by  changing  the 
geometry  parameters  try  to  change  frequency  range  of  interest.  First  idea  is  to  reduce  distance 
between  neighboring  screens.  In  Figures  4.19  and  4.20  results  for  several  screen  distances  are 
shown. 


SZmaxTE(0,0)ZmaxTE(0,0)dB 


Frequency  /  GHz 

Figure  4.19.  Sn  parameter  as  a  function  of  distance  between  neighboring  screens 
(distance  is  10mm,  9  mm,  8  mm  and  7  mm,  respectively) 
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SZminTM(0,0)ZmaxTE(0,0)dB 


distance  =  10 
distance  =  9 
distance  =  8 
distance  =  7 


Figure  4.20.  S21  parameter  as  a  function  of  distance  between  neighboring  screens  (coupling 
between  TE  input  and  TM  output  wave;  distance  is  10mm,  9  mm,  8mm  and  7  mm) 


With  smaller  distance  between  screens  realization  becomes  more  compact,  and  also  the 
working  frequency  is  shifted  to  the  lower  frequencies.  Second  idea  is  to  increase  dipole 
length  L.  With  larger  middle  metal  element  capacitance  between  edge  and  middle  element  is 
also  increased  and  shift  to  the  lower  frequencies  is  expected.  As  it  can  be  seen  in  Figures  4.21 
and  4.22,  the  operating  frequency  is  shifted  from  11  GHz  to  10  GHz,  and  the  bandwidth  is 
also  increased  from  1  GHz  to  almost  1.5  GHz. 


SZmaxTE(0,0)ZmaxTE(0,0)dB 


SZminTM(0,0)ZmaxTE(0,0)dB 


Frequency  l  GHz 

Figure  4.22.  S21  parameter  as  a  function  of  dipole  length  (coupling  between  TE  input  and  TM 

output  wave;  L= 8  mm,  8.3  mm  and  8.6  mm) 
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3rd  design 

In  order  to  increase  capacitance  some  changes  in  the  middle  element  geometry  could  be  made. 
Idea  is  to  increase  contact  surface  between  element  and  edge,  and  by  that  to  increase 
capacitance  between  load  and  surrounding  metal.  Therefore,  we  have  added  an  arch  to  the 
dipole.  By  changing  the  parameters  g,  r  and  B  we  can  affect  that  capacitance. 


bottom  of  dipole  are  introduced  to  increase  capacitance  (g=0.2, 13=45) 


SZmaxTE(0,0)ZmaxTE(0,0)dB 


SZrrinTM(0,0)ZmaxTE(0,0)dB 
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In  this  realization  there  is  a  shift  of  the  working  frequency  to  the  lower  frequencies  because  of 
increased  capacitance  between  the  dipole  element  and  surrounding  metal.  If  we  reduce 
distance  between  layers  from  10  mm  to  8  mm  we  get  Si i  which  is  unusable  in  practice. 


SZmaxTE(0,0)ZmaxTE(0,0)dB 


First  idea  how  to  improve  the  design  is  to  increase  parameter  g  (Fig.  4.23).  Increasing  the  g 
distance  between  arch  and  surounding  metal  is  decreased  and  capacitance  is  increased.  We 
can  expect  shift  to  the  lower  frequencies. 


SZmaxTE(0,0)ZmaxTE(0,0)dB 


=  0.1 
=  0.2 
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SZrninTM(0,0)ZmaxTE(0,0)dB 


Frequency  /  GHz 

Figure  4.28.  3rd  design:  S21  parameter  for  g=0. 1  mm,  0.2  mm  and  0.3  mm; 
coupling  between  TE  input  and  TM  output  wave  is  considered. 


By  increasing  angle  B  the  length  of  coupling  arc  between  middle  dipole  element  and  metal 
surface  is  increased.  Accordingly,  one  can  expect  increase  of  the  capacitance  and  shift  of  the 
working  frequency  to  the  lower  frequencies. 
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beta  =  20 
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Figure  4.29.  3rd  design:  Sn  parameter  for  B  =  20°,  30°,  40°  and  50° 


SZmnTM(0,0)ZmaxTE(0,C))dB 
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Figure  4.30.  3rd  design:  S21  parameter  for  B  =  20°,  30°,  40°  and  50°; 
coupling  between  TE  input  and  TM  output  wave  is  considered. 


4th  design 

The  basic  idea  of  the  4th  design  is  to  combine  two  orthogonal  loads  from  the  3rd  design.  By 
this  maybe  it  will  be  possible  to  increase  capacitance  and  get  lower  resonant  frequency,  and 
maybe  it  will  be  possible  to  obtain  polarizing  effect. 


Figure  4.31.  Double  cross  geometry 
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SZmaxTE(0,0)ZmaxTE(0,0)dB 


All  parameters  in  simulation  are  the  same  as  the  basic  parameters  in  third  realization.  As  we 
can  see  in  Figure  4.32  this  structure  is  useless.  Almost  all  energy  of  incident  wave  is  reflected, 
i.e.  it  is  not  possible  to  use  this  structure  as  a  polarization  rotator  or  a  polarizer. 
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CONCLUSIONS 
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Conclusions 


Metamaterial-based  cloaks  have  recently  been  proposed  to  prevent  scattering  of 
electromagnetic  waves,  i.e.  to  render  objects  invisible.  The  design  of  cloaks  is  usually  done 
using  transformation  electromagnetics,  i.e.  the  designed  structure  is  described  is  terms  of 
space-varying  permeability  and  permittivity  tensors.  However,  at  least  some  of  the  elements 
of  the  constitutive  tensors  should  have  close-to-zero  values  (either  relative  permeability  or 
permittivity).  Since  such  materials  do  not  exist  in  nature  one  should  realize  them  using 
metamaterial  approach  -  by  developing  an  appropriate  periodic  structure  with  period  smaller 
than  wavelength.  There  is  still  a  big  gap  between  the  wishes  -  desired  permeability  and 
permittivity  distribution  -  and  reality  -  the  real  electromagnetic  parameters  of  different 
metamaterial  structures.  The  purpose  of  this  project  is  to  make  a  step  in  filling  this  gap. 

Metamaterial  structures  are  extremely  complicated  to  design.  They  consist  of  a  large 
number  of  unit  cells  with  geometry  that  is  varying  along  the  considered  structure.  Due  to 
complexity  it  is  not  possible  to  design  them  without  good  strategy  (i.e.  we  cannot  simply  use 
optimization  possibilities  of  some  commercial  electromagnetic  solver).  Therefore,  we 
proposed  a  procedure  for  designing  metamaterial-based  devices.  The  procedure  starts  with  the 
ideal  design  obtained  using  some  physical  method  like  transformation  electromagnetics.  The 
ideal  design  contains  layers  of  homogeneous  anisotropic  materials  that  cannot  be  found  in 
nature.  Therefore,  we  need  to  search  for  metamaterial  realization  of  such  complex  multilayer 
structure.  In  order  to  do  that  in  an  efficient  way,  we  propose  an  intermediate  step  in  which  an 
equivalent  planar  anisotropic  multilayer  structure  is  defined  and  which  is  realized  using 
metamaterial  approach.  The  main  reason  for  introducing  an  intermediate  step  is  that  planar 
periodic  structures  can  be  rapidly  analyzed  (and  designed)  using  most  of  the  commercial 
electromagnetic  solvers  since  we  need  to  design  only  a  unit  cell.  After  designing  an 
equivalent  planar  structure,  we  need  to  go  back  and  map  (i.e.  to  curve)  the  structure  to  the 
desired  shape.  In  order  to  do  that  we  considered  an  efficient  approach  based  on  conformal 
mapping. 

We  have  also  investigated  methods  dealing  with  the  problem  of  how  to  electromagnetically 
characterize  metamaterial  structures.  The  first  idea  is  to  determine  effective  constitutive 
parameters  either  form  the  analysis  process  or  from  measurement  results.  However,  these 
methods  have  serious  drawbacks  if  the  considered  structure  is  resonant,  like  split-ring 
resonators.  In  that  case  it  is  not  possible  to  determine  effective  constitutive  parameters  around 
the  resonant  frequency  (unfortunately,  at  the  same  time  this  is  the  frequency  range  of  interest). 
Therefore,  we  propose  to  use  reflection  and  transmission  coefficients  for  characterization  of 
the  considered  structure,  which  (in  connection  with  the  optimization  routine)  represent  a 
powerful  tool  for  designing  metamaterial  structures. 
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We  have  also  investigated  the  method  that  can  simply  transform  the  auxiliary  planar 
structure  into  a  desired  shape.  It  is  shown  that  conformal  mapping  can  be  used  for  this 
purpose  resulting  in  accurately  determined  scattering  properties  of  the  final  structure.  We 
have  also  investigated  the  extension  of  the  Bloch  theory  for  multilayer  planar  periodic 
structure.  The  extension  to  cylindrical  structures  represents  an  accurate  method  for 
determining  properties  of  circular-cylindrical  multilayer  periodic  structures. 

Finally,  as  a  support  to  this  investigation  we  have  developed  several  algorithms  and 
programs.  First,  we  have  extended  the  previously  developed  G1DMULT  algorithm  to  include 
calculation  of  Green’s  functions  of  anisotropic  multilayer  structures  (the  G1DMULT 
algorithm  calculates  Green’s  functions  of  planar,  cylindrical  and  spherical  multilayer 
structures).  Furthermore,  we  have  developed  programs  for  analyzing  planar  and  cylindrical 
periodic  structures.  The  programs  are  based  on  the  Moment  Method  approach,  and  their  main 
purpose  is  fast  and  accurate  analysis  of  specific  periodic  structures. 

As  a  short  summary,  the  purpose  of  this  report  is  to  help  researchers  in  understanding 
where  the  critical  points  when  designing  metamaterial  structures  are,  and  to  propose  an 
efficient  and  successful  designing  procedure  of  metamaterial-based  devices.  We  sincerely 
hope  that  we  have  managed  to  give  contribution  in  that  direction. 
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APPENDIX 


A.l  DESCRIPTION  OF  G1DMULT  ALGORITHM 

The  algorithm  G1DMULT  calculates  spectral-domain  Green’s  functions  of  planar,  circular 
cylindrical  and  spherical  multilayer  structures  (see  [37]  for  details).  The  algorithm  will  be 
explained  for  the  planar  case.  The  3D  elements  are  replaced  by  equivalent  or  physical  current 
excitations  of  electric  and  magnetic  type,  which  then  are  Fourier  transformed.  In  the  planar 
case  we  use  the  two-dimensional  Fourier  transformation  in  the  x  and  y  directions,  defined  by 


f(kx ,  ky ,  z)  =  1 1  / (x,  y ,  z)e jxkx  e Jyky dxdy  (A.  1 ) 


By  this  Fourier  transformation,  the  3D  excitations  are  transformed  into  harmonic  current 
sheets.  If  the  source  is  infinitely  thin  in  z  direction,  we  get  one  discrete  current  sheet  per 
source,  otherwise  we  get  a  continuous  distribution  of  current  sheets  in  z-direction.  Each 
current  sheet  excites  two  plane  waves,  one  propagating  upwards  and  one  propagating 
downwards.  The  presence  of  the  multilayer  structure  will  cause  a  number  of  transmitted  and 

reflected  waves  in  each  layer,  with  variation  e~jk*x  and  e  Jk'  v  in  the  x-  and  y-directions. 
Notice  that  all  the  waves  have  the  same  field  variation  in  the  x-  and  y-directions  for  all  layers, 
which  is  a  consequence  of  the  boundary  conditions.  Therefore,  only  the  field  variation  in  the 
direction  perpendicular  to  the  boundaries  is  unknown,  so  we  have  a  harmonic  one¬ 
dimensional  (ID)  field  problem.  In  this  way  the  spectral  domain  problem  is  interpreted  as  a 
ID  spatial  domain  problem  consisting  of  the  ID  multilayer  structure  and  harmonic  ID 
sources  in  the  form  of  current  sheets. 

The  algorithm  G1DMULT  is  based  on  subdividing  the  spatial  harmonic  problem  into  one 
equivalent  problem  per  layer,  where  the  field  in  each  region  is  formulated  as  the  field  radiated 
by  equivalent  currents  at  the  layer  boundaries.  For  example,  the  E-field  in  the  layer  j  is 
expressed  as 


E,= 


Ghomj  .  horn  t  . 

FI  J  /—I  +  vJ  RJ  J  v  + 


homi 


'  EJ 


’7-1 


EJ 


+ 


hom  t  exci 


+  G“mJ 


+ 


Gh°rm\ /rex< 
EM 


(A.2) 


where  J;  and  M .  are  equivalent  electric  and  magnetic  current  sheets  at  boundary  j,  J "" 
and  M.exci  are  excitation  electric  and  magnetic  currents  in  layer  j  (if  any),  and  Ghomis  the 
Green's  function  of  the  homogeneous  problem.  By  using  J .  =  ±  z,  x  FL  and 
M  ■  =  +  |xE.  eq.  (A.2)  can  be  expressed  in  terms  of  the  unknown  tangential  E-  and  H- 
fields  E  .  and  H  at  the  boundary  j  between  layers  j  and  j+1  and  known  excitation  currents. 
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The  boundary  conditions  that  the  tangential  E-  and  H-  fields  are  continuous  at  the  layer 
boundaries  give  4  linear  equations  per  boundary.  The  tangential  E-  and  H-fields  are  evaluated 
by  solving  the  system  of  4(Niayer-l)  equations  with  4(Niayer-l)  unknowns,  where  Niayer  is  the 
number  of  layers  in  the  multilayer  structure.  After  they  have  been  determined,  the  total  E-  and 
H-fields  at  any  desired  z-location  can  be  found  by  using  the  equivalent  problem  formulation. 

The  core  problem  in  the  formulation  is  to  calculate  the  E-  and  H-  fields  due  to  a  harmonic 
current  sheet  in  a  homogeneous  region.  The  two  subroutines  G1DJ  and  G1DM  which 
calculate  these  fields  (G1DJ  for  electric  source  and  G1DM  for  magnetic  source)  are  the  only 
part  of  the  routine  which  depend  on  whether  we  consider  a  planar,  circular  cylindrical  or 
spherical  geometry,  and  they  represent  the  only  difference  between  the  three  versions  of  the 
GIDMULTroutine.  In  the  planar  case,  the  electromagnetic  field  excited  by  a  current  sheet  is 
equal: 


E  = 


/7H  = 


2k 


2k 


1 - 1 
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(A.2) 


(A.3) 


where  k+ =(kxx  +  kyy  +  kzz)/ k  and  k  =(kxx  +  k  y  —  kzz)/  k .  In  the  cylindrical  case,  the 
electromagnetic  field  excited  by  a  current  tube  is  equal : 


(a)  from  z-directed  current  tube 


P^Ps 
P^Ps 

Hz(p,n,kz)  =  0 

(b)  from  <j>  -directed  current  tube 


Ez(p,n,kt)  =  ~^P,  J-Sn.k,) 


»?’(*„ A)  JAKP) 


K(KP.)  H"\ko) 


(A.4) 

(A.5) 


E(nnk)-  35*1  j(nk).\H"  JJk”p)  p-p- 

ESp.n.k,)  -  2  k  H?{k  p>_pi 


(A.6) 


H  (p  n  k  )  -  -—k  p  Jink)-  K’Aa>  '.<V>  P*P. 
-  P  2  UsKp.)  H?(kpP)  pip, 


(A.7) 


The  other  field  components  are  calculated  using  the  following  expressions: 
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EP  = 
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More  details  about  the  G1DMULT  algorithm  can  be  found  in  [37]. 


(A.  8.  a) 


(A.8.b) 


(A.8.c) 


(A.8.d) 


Distribution  A:  Approved  for  public  release;  distribution  is  unlimited. 


Constitutive  Parameters  of  Metamaterial  Structures  Used  for  Invisible  Cloak  Realization  69 


A.2  GREEN’S  FUNCTIONS  OF  ANISOTROPIC  PLANAR  STRUCTURE 


Let  E.  H  :  M3  — >  C3  represents  electric  and  magnetic  field  respectively.  Components  of 

mentioned  fields  will  be  annotated  E  —  x  •  E  where  hat  over  x  means  unit  vector  in  x 

direction.  In  the  case  of  partial  derivation  over  some  variable,  notation  will  be  slightly 
changed.  Component  mark  in  that  case  will  be  in  superscript,  and  derivative  mark  will  be  in 

subscript  Ex  =  d  E  =d  X  •  E  .  Fourier  transform  will  be  annotated  with  hat  over 

y  y  x  y 

function.  In  this  appendix  next  definition  of  Fourier  transform  will  be  used: 


E  kx,ky,z  —  J  E  x,y,z  e  ^ yVdxdy 

M2 

E  x,y,z  = — E  ejkxXejkyVdkxdky 

2n  M2 


(A.9) 


A2.1  Case  1:  Source  with  direction  orthogonal  to  the  axis  of  anisotropy 

In  the  first  case  current  is  orthogonal  to  the  axis  of  anisotropy.  Permeability  and  permittivity 
will  be  diagonal  matrices  with  same  entry  at  xx  and  yv  places,  and  different  at  zz.  Because  the 
axis  of  anisotropy  is  in  the  z  direction,  all  other  field  components  can  be  described  over  z 
components  of  electromagnetic  field.  Therefore,  the  idea  is  to  find  differential  equations  for  z 
component  of  electric  and  magnetic  field.  Permeability  and  permittivity  looks  like: 


dx 
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(A.  10) 


In  considered  situation  excitation  current  will  have  only  x  and  y  components: 


J  = 


axS(x,  y,  z) 
x , 

0 


ay6{x,y,z) 


From  Maxwell  equations: 

V  x  H  =  iuj£QeE 


(A.ll) 
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we  get  set  of  equations  for  electric  field: 

V  x  /i-1V  x  E  —  kfisE  =  0,  =  u2l~i0£0  (A.  12) 

Because  axis  of  anisotropy  is  z  axis,  we  only  need  equation  for  z  component  (third  equation): 

Ex  -  Ez  -Ez  +  Ey  -  khi  e  E  =0  (A.13) 

zx  xx  yy  1  zy  0  r^x  z  z  v  7 

Using  the  Gauss  law  for  electrical  field  div  £q£E  —  0  we  get  equation  which  connects  all 
three  components  of  electric  field: 

+  CyEl  +  £zEl  =  0  <A'14> 

From  that  equation  we  can  express: 

Ex+Ey  =  -^EZ  (A.  15) 

x  1  y  z  v  7 

£ 

^x 

From  equations  (1)  and  (2)  we  can  express  equation  for  z  component  of  electric  field: 

Ez  +  Ez  +£ez  +  k2u  £  E  =0  (A.  16) 

xx  1  yy  '  zz  1  O^x  z  z  v  7 

ex 

Using  the  duality  theorem  we  can  get  equation  for  z  component  of  magnetic  field: 

H‘xx  +  H;v+^~  HI  +  %ex»zHz  =  0  (A.17) 

rx 

This  equation  is  partial  differential  equation  in  all  three  space  variables.  That  equation  is  not 
standard  Helmholtz  equation  because  it  contains  different  constants  multiplying  second 
derivative  of  field  component.  The  idea  is  to  use  Fourier  transformation  over  x  and  y  variables 
and  transform  equations  (3)  and  (4)  into  ordinary  differential  equations  over  variable  z: 

~k2A  -  eyEz  +£-^Ezzz+  kfax£zEz  =  0  (A- 18) 

£ 

°x 

This  is  equivalent  to  Helmholtz  equation  for  the  z  component  of  electrical  field: 

El,  +  klEz  =  0  kl=  k&A  --01-  (A.19) 

^  Z 

From  the  duality  theorem  we  get: 

K+k\Hz=  0  kl=^xex-^P\  (A.20) 

^Z 

where  02  =  £  +  £ . 

Two  independent  solutions  of  the  Helmholtz  differential  equation  are  combinations  of 
complex  exponentials.  From  the  Silver-Muller  radiation  condition  we  know  that  solutions 
have  this  form: 
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Ae  jkE\z  >  0 
BejkEZ  ,z  <  0 


Ce  jk*z,z  >  0 
Dejknz  ,z  <  0 


(A.21) 


All  other  components  of  electric  and  magnetic  field  will  be  expressed  over  z  components. 
From  transformed  Ampere  and  Faraday  law  we  can  express  other  field  components  in  the 
form  of  Helmholtz  equation  with  nonzero  right  side: 


kWA  +  El 
+  ei 
+  hi 
*oVA  +  Hi 


“fW’JtyH,  +  AE* 
-“WxKH,  +  AEz 
Ah*  -  “o eAEz 

AA  +  4 


(A.22) 


From  the  general  form  of  solution  for  z  components  of  electric  and  magnetic  field,  and  from 
above  equations  we  can  find  expressions  for  all  other  components: 
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Boundary  conditions  for  our  situations  are: 


K~Hp=ay 


H+ 
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Hv  = 
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E+  -  E~  =0 

x,y  x,y 


(A.23) 


(A.24) 


By  fulfilling  the  boundary  conditions  we  get  coefficients  for  the  z  components  of  electric  and 
magnetic  field: 


C  =  D  = 


— — —  k  a  —  k  a 
2 kH/iz  x  y  y  x 
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A2.2  Case  2:  Source  with  direction  parallel  to  the  axis  of  anisotropy 

In  the  second  case  the  source  current  is  parallel  with  the  axis  of  anisotropy.  Permeability  and 
permittivity  will  be  diagonal  matrices  with  same  entry  at  xx  and  zz  place,  and  different  at  yy. 
Since  the  axis  of  anisotropy  is  in  the  y  direction,  all  other  field  components  can  be  described 
over  y  field  components.  Permeability  and  permittivity  tensors  look  like: 
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In  the  described  situation  current  will  have  only  y  component  of  the  source: 


J  = 


0 

ayS(x,y,z) 

0 


(A.27) 


(A.28) 


Since  the  source  current  is  in  the  y  direction,  magnetic  field  in  y  direction  is  equal  to  zero. 
Because  of  that  all  other  components  of  electric  and  magnetic  field  could  be  expressed  over  y 
component  of  electric  field.  From  set  of  equations  (A.22)  we  use  only  the  second  equation 
(equation  for  y  component): 

Ez  -Ex  -Ey  -Ey  -kill  e  E  =0  (A.29) 

yz  yx  xx  zz  0  >xy  y  v  ’ 

Gauss  Law  gives  us  equation  which  connects  all  electric  field  components,  and  with  that 
equation  we  find  final  equation  for  y  component  of  electric  field: 

EL+-  El  +  El  +  kW,E,  = 0  <A'30> 

£x 

In  spectral  domain  we  get  Helmholtz  equation: 

Eyy  +  k2yEy  =0  k2y  =  %lx£x  -^P2  (A.31) 

£y 

where:  /32  =  k2  +  k2 . 


Solution  of  that  equation  have  a  form: 
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(A.32) 
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All  other  components  of  electric  and  magnetic  field  will  be  expressed  over  y  components. 
From  transformed  Ampere  and  Faraday  law  we  can  express  other  field  components  in  the 
form  of  Helmholtz  equation  with  nonzero  right  side: 


klu  e  E  +EX  =  jk  Ey 

O^V x  x  1  yy  J  x  y 


kin  e  E  +EZ  =  jk  Ey 

0  r“x  x  z  1  yy  J  z  y 


ki +  Hyy  =  u£0£xKE; 


(A.33) 


k0^£xHz  +  Hyy  =  -U£0£xKEy 


From  the  general  form  of  solution  for  y  components  of  electric  and  magnetic  field,  and  from 
above  connection  equations  we  can  find  expressions  for  all  other  components: 
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H  = 


^fof yK  £ 
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In  the  considered  situation  there  are  only  two  undetermined  coefficients.  They  can  be 
determined  using  definition  of  equivalent  magnetic  current: 


Meq  =  E"1— V*  J 


(A.35) 


UJ£. 


0 


In  our  case  electric  current  source  has  only  y  component,  which  is  transformed  with  Fourier 
transform  into: 


Kq  =  £ 


We  will  use  only  two  boundary  conditions: 
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(A.36) 


h;  =  o 


E+  -  E~  =  -M. 


(A.37) 


From  those  two  boundary  conditions  and  from  definition  of  equivalent  magnetic  current  we 
find  undetermined  coefficients: 


A  =  B  = 


Bl 


(A.38) 
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A3.  GREEN  S  FUNCTIONS  OF  ANISOTROPIC  CYLINDRICAL  STRUCTURES 

In  the  considered  case  permeability  and  permittivity  tensors  will  be  diagonal  matrices  with  the 
same  entry  at  pp  and  <j)<P  place,  and  with  different  entry  atzz  .  Because  the  axis  of  anisotropy 
is  in  the  2  direction,  all  other  field  components  can  be  described  over  z  components  of 
electromagnetic  field.  Permeability  and  permittivity  tensors  looks  like: 
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In  the  considered  situation  the  source  current  will  have  only  2  component: 


Starting  equations  are: 


J  =  zL 


S(P  -  d) 

2i vp 


V  x  n  1V  x  E  —  k^eE  =  0,  k^  =  cj2/j,0£0 

div(eE )  =  0 


(A.39) 


(A.40) 


(A.41) 


From  the  definition  of  nabla  operator  in  the  cylindrical  coordinate  system  we  get  equation  for 
2  component  of  electric  field: 

-  - 1 rr  -  = 0  <A42) 

£p  dzz 

A  represents  Laplace  operator  over  p  and  </>  variables.  For  the  reason  of  easier  computation 
some  abbreviations  will  be  used.  Laplace  operator  will  be  decomposed  into  two  operators. 

-A  ,  =  A  +—A.  (A.43) 

P<P  P  <P  v  ' 

The  idea  is  to  use  the  Fourier  transformation  over  two  variables  to  transform  the  considered 
partial  differential  equation  into  an  ordinary  differential  equation.  The  differential  operator  A. 

will  be  transformed  into  multiplication  with  a  constant.  Our  cylinder  is  infinitely  long  in  the  z 
direction  and  2ir  periodic  in  <f>  variable.  Because  of  that  the  ordinary  Fourier  transform  in  z 
direction  will  be  used: 


“  K 

U  X 


(A.44) 


In  </>  variable  the  ordinary  Fourier  series  will  be  applied: 
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a. 


2tt 

= -kJf(x)e~’' 


Jnxdx 


(A.45) 


f(x) 


V27T 
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ri  — — oo 


ane 


jnx 


Using  this  two  transformations  we  get  an  ordinary  differential  equation  in  p  variable: 


ld_ 

pdp 


/  ^  \ 

dEn 

(  2  ^ 
x  9  n 

p  z 

+ 

x 2  -  — 

dp 

\  '  J 

p  J 

E™  =  0  with  A  =  k*ii  ez  — -k 


er 


(A.46) 


This  is  the  Bessel  equation  of  order  n.  Solution  of  this  equation  is  coefficient  E™  in  Fourier 
expansion  of  E  field  component.  Physical  conditions  that  we  have  to  satisfy  are: 


lim  E™(p,kz)  <  oo  lim  E"(p,  kz)  =  0  (A.47) 

p— >  oo  p  — >0 


Because  of  condition  (A-47)  solution  will  be  represents  with  Bessel  function  inside  cylinder 
and  with  outgoing  Hankel  function  outside  cylinder  with  radius  d. 


AnJJXP)  iP<d 
BnH(P(\p),p  >  d 


p-K 


CnJnM  ,P<d 

DnH^(Xp),p  >  d 


(A.48) 


From  Maxwell  equations,  from  two-dimensional  Fourier  transformation,  and  by  using 
notation  £2  —  kip  e  —  k2  we  get: 

0 r  p  p  Z  ° 


1  fln  ,  3\  dp 

f!  p  z  0  dp 

wyv  9H^  nkz  1 

Q  dp  Q  p  z 

dp  Q  p  z 

nK  1  pn  dE; 

Q  p  z  Q,  dp 


(A.49) 


Boundary  conditions  in  the  considered  situation  are: 


Distribution  A:  Approved  for  public  release;  distribution  is  unlimited. 


Constitutive  Parameters  of  Metamaterial  Structures  Used  for  Invisible  Cloak  Realization  76 
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(A.50) 


By  using  the  boundary  conditions  in  cylindrical  coordinate  system,  Wronskian  for  Bessel 
functions,  definition  of  other  field  components  over  z  components  and  solution  for  the  z 
component  of  electric  and  magnetic  field  we  get  expressions  for  coefficients  in  the  considered 
expansion: 


K  pA 
h:  pa 


Jbrd_jn 

2  Ld£n£_  Z 


=  l- 


"Op 

.  Q,7rd2X 


2  nk„ 


H^(Xd)Jn{Xp),p  <  d 
Jn(\d)H®(Xp),p  >  d 
H^(Xd)'Jn(Xp),p  <  d 
JJXd)'H^(Xp),p  >  d 


(A.51) 


By  using  the  definition  of  inverse  Fourier  transform  over  p  and  6  variables  we  get  the  final 
solution: 


EzAy,z) 
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(A.52) 
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A.4  DESCRIPTION  OF  MOMENT  METHOD  PROGRAM  FOR  ANALYZING 
PLANAR  AND  CIRCULAR-CYLINDRICAL  PERIODIC  STRUCTURES 


Periodic  strips  can  be  accurately  analyzed  by  expanding  the  unknown  currents  on  metal 
patterns  in  basis  functions,  and  by  using  MoM  to  numerically  determine  the  amplitudes  of  the 
basis  functions.  We  formulate  an  integral  equation  by  stating  that  the  tangential  electric  field 
is  zero  at  the  metal  surface,  i.e.  hx  [e'"c  +EiCO']=  0.  Since  the  structure  is  periodic  the 
electromagnetic  field  excited  by  the  currents  is  in  the  form  of  the  Floquet  modes,  and  the 
currents  on  different  metal  patterns  are  identical  except  for  a  linear  phase  difference  equal  to 
the  considered  harmonic  variation  of  the  incident  field. 


The  problem  is  solved  in  the  spectral  domain.  In  the  planar  case,  we  use  a  two-dimensional 
Fourier  transformation  in  the  directions  in  which  the  structure  is  periodic  (i.e.  x-  and  y- 
directions)  and  in  the  cylindrical  case  we  use  a  Fourier  transformation  in  a  axial  z  direction 
and  a  Fourier  series  in  a  circumferential  (j)  direction.  The  scattered  field  Escat  can  be  expressed 

in  terms  of  the  spectral  domain  dyadic  Green's  function  G  and  the  Fourier  transform  of  the 
current  on  one  metal  pattern,  according  to 


Wcat(un,u0,u,)  = 


P  P. 


X  X  gk. i o-iK 


1  \  „-JkoUo  ~KUP 


Xp)e 


m=- oo  /=— oo 


K 


2m  n 


+  k : 


21k 


+  k 


inc 

P 


(A.53) 

(A.54) 


where  k“cand  k"'c  are  the  propagation  constants  of  the  incident  wave  in  two  orthogonal 

directions  across  the  metal  pattern  (i.e.  o  and  p  directions),  P0  and  Pp  are  the  periodicity  of  the 
unit  cell  in  o  and  p  directions,  and  '  denotes  the  source  coordinates.  The  spectral  domain 

Green's  function  Gis  calculated  by  the  G1DMULT  routine.  The  integral  equation  is 
transformed  into  a  matrix  equation  by  using  the  MoM  where  the  test  functions  are  the  same  as 
the  basis  functions  (Galerkin's  method),  and  the  inner  product  is  taken  across  the  central  unit 
cell.  The  expression  for  the  elements  of  the  MoM  matrix  is  derived  from  eq.  (A.53),  and  the 
expression  for  the  elements  of  the  voltage  vector  is  obtained  by  integrating  the  product  of  the 
test  function  and  the  incident  field  component  parallel  to  the  test  function.  The  elements  of  the 
impedance  matrix  and  the  excitation  vector  are 


00  00 


Z,  = 


P  P 

1  1  ,,  m= 


X  X^-k 


o  p 


m=-oo  /=— oo 


V,  =  J(-CV*D '  E (u\ ,  u\  =  0,  u\  =  0) 


(A.55) 

(A.56) 


where  conveniently  we  choose  that  the  w'0and  u  p  coordinates  of  the  center  of  the  metal 
pattern  are  zero.  In  the  planar  case,  the  quantities  in  eqs.  (A.53)-(A.56)  are: 


un  =  z>  u0=x>  uP  =  y 

In  the  cylindrical  case,  the  incident  plane  wave  is  expanded  as 


(A.57) 
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E‘znc  =  E0  sin  0inc  cos  ainc  ■  £  j"Jn  ( k0p sin  (A.58) 

n=- oo 

where  J„  is  the  nth-order  Bessel  function  of  the  first  kind,  the  polarization  angle  amc  is  the 
angle  the  incident  electric  field  makes  with  the  plane  of  incidence,  and  &nc  and  <j)mc  are  angles 
of  incidence  measured,  respectively,  from  the  z-axis  and  x-axis  in  the  xy  plane.  The  scatterer 
has  the  following  symmetry:  by  rotating  the  scatterer  by  2 nj N rji  where  Nt/)  is  the  number  of 

unit  cells  in  (j)  direction,  we  get  the  same  structure.  Thus  we  expand  the  fields  in  Floquet 
modes  in  <f>  direction  where  for  every  component  of  the  incident  field  in  eq.  (A.58)  the  strips 
excite  fields  which  have  the  same  phase  variation  between  the  centers  of  the  strips  as  the 
incident  cylindrical  wave.  In  this  case  the  scattered  field  due  to  the  incident  wave  with 
exp (—jn</>)  variation  is  also  given  by  eq.  (A.53)  with 

un  =  A  K  =  (j),  up=z  (A.59) 

and 

k‘0nc=n,  k™c  =  -k[nc  =  -k0  cos  9inc .  (A.60) 
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